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ABSTRACT 


We present Atacama Large Millimeter/submillimeter Array (ALMA) sub-kiloparsec- to kiloparsec-scale resolution observations of 
the [Cm], CO (9-8), and OH* (1,—01) lines along with their dust continuum emission toward the far-infrared (FIR) luminous quasar 
SDSS J231038.88+185519.7 at z = 6.0031, to study the interstellar medium distribution, the gas kinematics, and the quasar-host 
system dynamics. We decompose the intensity maps of the [Cm] and CO (9-8) lines and the dust continuum with two-dimensional 
elliptical Sérsic models. The [Cu] brightness follows a flat distribution with a Sérsic index of 0.59. The CO (9-8) line and the dust 
continuum can be fit with an unresolved nuclear component and an extended Sérsic component with a Sérsic index of ~1, which may 
correspond to the emission from an active galactic nucleus dusty molecular torus and a quasar host galaxy, respectively. The different 
[C rr] spatial distribution may be due to the effect of the high dust opacity, which increases the FIR background radiation on the [Cm] 
line, especially in the galaxy center, significantly suppressing the [C rr] emission profile. The dust temperature drops with distance 
from the center. The effective radius of the dust continuum is smaller than that of the line emission and the dust mass surface density, 
but is consistent with that of the star formation rate surface density. This may indicate that the dust emission is a less robust tracer 
of the dust and gas distribution but is a decent tracer of the obscured star formation activity. The OH* (1;—0,) line shows a P-Cygni 
profile with an absorption at ~-400 km/s, which may indicate an outflow with a neutral gas mass of (6.2 + 1.2) x 105 Me along the 
line of sight. We employed a three-dimensional tilted ring model to fit the [Cm] and CO (9-8) data cubes. The two lines are both 
rotation dominated and trace identical disk geometries and gas motions. This suggest that the [Cm] and CO (9-8) gas are coplanar 
and corotating in this quasar host galaxy. The consistent circular velocities measured with [Cu] and CO (9-8) lines indicate that 
these two lines trace a similar gravitational potential. We decompose the circular rotation curve measured from the kinematic model 
fit to the [C n] line into four matter components (black hole, stars, gas, and dark matter). The quasar-starburst system is dominated 
by baryonic matter inside the central few kiloparsecs. We constrain the black hole mass to be 2.97193) x 10° Mo; this is the first 
time that the dynamical mass of a black hole has been measured at z ~ 6. This mass is consistent with that determined using the 
scaling relations from quasar emission lines. A massive stellar component (on the order of 10? Mo) may have already existed when 
the Universe was only ~0.93 Gyr old. The relations between the black hole mass and the baryonic mass of this quasar indicate that 
the central supermassive black hole may have formed before its host galaxy. 
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1. Introduction 


Almost 400 quasars at z > 5.7 have been discovered in the 
past ~20 years, mostly from optical wide-field multiband sur- 
veys (e.g., Fan et al. 2000; Matsuoka et al. 2022). These quasars 
provide a unique opportunity to study a number of key issues, for 
example the formation of young luminous quasars, the evolving 
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impact of the central black holes on the host galaxies, and the 
typical interstellar medium (ISM) conditions in the quasar host 
galaxies, during the epoch at which the intergalactic medium was 
being reionized by the first luminous sources. 


Toward some of these z ~ 6 quasars, bright millimeter dust 
continuum emission has been detected at millijansky or sub- 
millijansky levels (e.g., Bertoldi et al. 2003a; Wang et al. 2007, 
2008, 2011a; Omont et al. 2013; Venemans et al. 2018; Li et al. 
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2020c) using the IRAM facilities, the James Clerk Maxwell Tele- 
scope (JCMT), and the Atacama Large Millimeter/submillimeter 
Array (ALMA), indicating far-infrared (FIR) luminosities of 
10!! to 10!3 Lo and star formation at rates of hundreds to thou- 
sands Mo/yr in these young quasar host galaxies. A few of these 
z ^ 6 FIR luminous quasars have been detected in multi-J CO 
transitions from J = 2-1 to 17—16 (e.g., Bertoldi et al. 2003b; 
Walter et al. 2003, 2004; Carilli et al. 2007; Riechers et al. 2009; 
Wang et al. 2010, 2011a,b; Gallerani et al. 2014; Shao et al. 
2019; Li et al. 2020b; Pensabene et al. 2021) using the IRAM 
facilities and ALMA, indicating abundant gas reservoirs on the 
order of 101? Mo and highly excited molecular gas in the quasar 
host galaxies. Most of these high-z FIR-luminous quasars have 
been detected in [C n] 158 um fine-structure line emission down 
to arcsecond- and subarcsecond-scale resolution using the IRAM 
facilities and ALMA (e.g., Maiolino et al. 2005; Walter et al. 
2009; Wang et al. 2013; Shao et al. 2017; Decarli et al. 2018; 
Izumi et al. 2018; Venemans et al. 2020; Yue et al. 2021; Walter 
et al. 2022; Meyer et al. 2022), suggesting that high-luminosity 
quasar host galaxies have lower dynamical masses than local 
galaxies with similar black hole masses, although the ratios of 
the black hole mass to the host galaxy dynamical mass of most 
of the low-luminosity quasars are consistent with the local value. 
In addition, high-resolution CO, [C u], and ultraviolet (UV) con- 
yy tinuum observations reveal close companions for some z ~ 6 
^ quasars (e.g., Wang et al. 2011b, 2019; McGreer et al. 2014; 


» Decarli et al. 2017; Miller et al. 2020; Izumi et al. 2021; Pens- 


» abene et al. 2021). These studies suggest an overall scenario of 
. the coevolution of quasars and their host galaxies: the central 
supermassive black holes (SMBHs) may grow rapidly through 
major galaxy mergers and the accretion of large amounts of gas, 


J which may trigger high rates of star formation in the quasar host 
YJ galaxies and thus influence the relationship between the black 


hole mass and the dynamical mass of these quasar-host systems. 

However, only a few of these studies are based on high- 
| resolution (i.e., sub-kiloparsec scale) and high sensitivity ob- 
. servations in the millimeter (e.g., Yue et al. 2021; Walter et al. 
2022). The quasar host galaxies at z ~ 6 are at most a few beam 
sizes across. Thus, the inferred ISM properties and the dynami- 


© cal information of the quasar-host systems are mostly based on 


unresolved or slightly resolved spatially integrated flux stacking 
along the velocity or frequency axes (i.e., the intensity maps) and 
on the distribution of spatially integrated flux for each velocity 
) interval as a function of apparent radial Doppler velocity (i.e., the 
integrated spectra). High-resolution (i.e., a few hundred parsecs) 
ALMA observations of the ISM in these quasar host galaxies at 
z ~ 6 can provide insight into these issues, for example the spa- 
tial distribution of different ISM tracers in addition to the star 
formation activity, the gas kinematics, and the overall dynamics 
probed by different ISM tracers in the SMBH-host system at the 
earliest epochs. 

Measuring the sizes and morphologies of galaxies in the 
early Universe is critical for understanding the initial stage of 
galaxy formation and evolution. The long-wavelength FIR dust 
continuum traces regions of young and compact star formation 
that are severely obscured by dust at optical wavelengths. In both 
observed galaxies (e.g., Riechers et al. 2013, 2014; Ikarashi et al. 
2015; Hodge et al. 2016; Wang et al. 2019; Tadaki et al. 2020) 
and simulated galaxies (e.g., Cochrane et al. 2019; Popping et al. 
2022), the distribution of the dust continuum emission generally 
is more compact than both the cold gas and the dust mass, but 
is more extended (more compact) than the stellar component at 
z € 3 (z > 4). The [C u] fine-structure transition at 157.74 um is 
one of the primary coolants of the star-forming ISM in galaxies 
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(e.g., Stacey et al. 1991), and it traces both the neutral atomic and 
ionized gas. The different rotational transitions of CO trace cold 
and warm gas that will form stars in the future. Their sizes and 
morphologies are strong probes of ISM properties. Thus, com- 
paring the size and morphology between the dust emission and 
the [C n] and CO lines will provide insights into the evolutionary 
process of the ISM in which star formation takes place. 

The kinematics of the ISM can yield important information 
on the dynamical structure of the quasar host, as well as the dis- 
tribution and fraction of baryonic and dark matter in the system 
(e.g., van de Hulst et al. 1957; Rubin 1983; de Blok et al. 2008). 
The rotation velocities of different phases of the ISM frequently 
show a complex picture and differ from source to source. The ro- 
tation velocities measured from the molecular components (i.e., 
CO) and atomic components (i.e., H1) are generally in agree- 
ment (e.g., Wong et al. 2004; Frank et al. 2016). Übler et al. 
(2018) find that the kinematics of CO and H a are in good agree- 
ment in a z — 1.4 star-forming galaxy. However, Simon et al. 
(2005) find that the H a line in NGC 4605 shows systematically 
slower rotation than the CO. de Blok et al. (2016) find that the 
velocity of the [C n] line is systematically larger than that of the 
CO or H t lines in a few nearby galaxies, which they attribute to 
systematics in the data reduction and the low velocity resolution 
of the [Cm] data. It is still not clear why discrepancies between 
different kinematic tracers are observed in some galaxies. 

In the nearby Universe, the stellar bulge can be observed di- 
rectly through optical and near-infrared (NIR) imaging, and its 
dynamics can be probed via investigations of the rotation curve 
(e.g., Begeman et al. 1991; de Blok & Bosma 2002; Sofue et al. 
2009; Gao et al. 2018). At high redshifts, the stellar bulge may 
already exist long before the peak of cosmic star formation, as 
demonstrated by the [Cu] gas kinematics of 4 < z < 5 dusty 
star-forming galaxies and quasars (e.g., Rizzo et al. 2020, 2021; 
Tsukui & Iguchi 2021). However, at higher redshifts, into the 
reionization epoch, it remains unknown whether the stellar bulge 
is already present. 

In this paper we report on ALMA sub-kiloparsec- to 
kiloparsec-resolution observations of the [Cu], CO (9-8), and 
OH? (1,—01) lines and their underlying dust continuum toward 
the FIR-luminous quasar SDSS J231038.88+185519.7 (here- 
after J2310+1855) at z = 6.0031, to study the ISM distribution, 
the gas kinematics, and the quasar-host system dynamics. This 
quasar is first reported by Wang et al. (2013) with 250 GHz dust 
continuum and CO (6—5) line observations using the IRAM fa- 
cilities and [C n] line observations using the ALMA. Jiang et al. 
(2016) are credited with its discovery; they discovered this opti- 
cally bright quasar with m,4594 = 19.30 mag in the Sloan Digital 
Sky Survey (SDSS) imaging data. The black hole mass has been 
measured to be (4.17 + 1.02) x 10? Mo and (3.92 + 0.48) x 10? Mo 
from the GEMINI/GNIRS spectra of Mg and C v lines, re- 
spectively (Jiang et al. 2016). A smaller black hole mass of 
(1.8 + 0.5) x 10? M, based on the C iv emission line detected 
in the X-SHOOTER/VLT spectrum is reported by Feruglio et al. 
(2018). All of these estimates of the black hole mass are based on 
the local scaling relations (e.g., Shen 2013). However, it is still 
under debate if the local relationship is suitable at high redshift. 

J2310+1855 has a well-measured dust content and a detailed 
rest-frame NIR-to-FIR spectral energy distribution (SED) ob- 
served by SDSS, WISE, Herschel, IRAM, and ALMA. (Shao 
et al. 2019). These observations reveal a dust temperature of ^40 
K and a star formation rate (SFR) of 2000 Mo/yr under the op- 
tically thin approximation. However, the dust and star formation 
distribution within the host are unknown, and the origins of dif- 
ferent dust components have not yet been identified. The multi-J 
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CO transitions (from J = 2-1 to 13-12) have been detected using 
the Karl G. Jansky Very Large Array (VLA), the IRAM facili- 
ties, and ALMA (Wang et al. 2013; Feruglio et al. 2018; Shao 
et al. 2019; Carniani et al. 2019; Li et al. 2020b; Riechers et al. 
in preparation). A detailed CO spectral line energy distribution 
analysis reveals that, in addition to the far-UV radiation from 
young and massive stars, another gas heating mechanism (e.g., 
X-ray radiation and/or shocks) may be needed to explain the ob- 
served CO luminosities (Carniani et al. 2019; Li et al. 2020b). 
The low excitation water para-H5O (202-111) and para-H5O (325- 
313), OH* (11—01) lines, and lines from ionized gas such as [O 1] 
146 um, [O m] 88 um, and [N u] 122 um have been detected with 
ALMA toward J2310--1855, and a solar-level metallicity is pro- 
posed based on the [O m]/[N u] ratio (Hashimoto et al. 2019b; 
Li et al. 2020a; Tripodi et al. 2022). We should note that [N u] 
122 um is an upper [N u] line (i.e., the transition of 3p53p4), 
so there are potential excitation effects. In addition, [O m] 88 um 
and [N ir] 122 um lines have significantly different ionization po- 
tentials (^14 versus ~35 eV), so they do not necessarily trace the 
same parts of the H n regions. The bright [C n] line has been de- 
tected with ALMA at low angular resolution (> 077; Wang et al. 
2013; Feruglio et al. 2018). A velocity gradient is obvious, which 
may indicate that the disk is dominated by rotating gas. The ~4 
kpc resolution [C rr] data reveal a dynamical mass of 9.6 x 10!° 
Me with an approximate estimate of the inclination angle (46°, 
determined from the ratio between the minor and major axis), 
suggesting a Mgy/Mpulge value that is higher than the local value 
(Wang et al. 2013). Tripodi et al. (2022) used ~1 kpc resolution 
[C u] data to find a dynamical mass of 5.2 x 10!° Mo within a 1.7 
kpc region, based on a kinematic modeling on the data cube with 
a rotating disk inclination angle of 25°. However, the limited spa- 
tial resolution introduces large uncertainties in the determination 
of the gas kinematics, making it difficult to perform a dynamical 
decomposition of the quasar-host system. In this paper we use 
ALMA high-resolution observations of the [C i] and CO (9-8) 
lines (~0.6 and ~1 kpc, respectively) and their underlying dust 
continuum in this quasar to explore its dynamics in detail. 

The outline of this paper is as follows. In Sect. 2 we describe 
our ALMA observations and the data reduction. In Sect. 3 we 
; present the measurements of the [C rn], CO (9-8), and OH" (1,- 
01) lines and the underlying dust continuum. In Sect. 4 we fit 
the line and continuum intensity maps to 2D Sérsic functions, 
constrain the dust properties, apply a 3D tilted ring model to 
the [Cu] and CO (9-8) data cubes, and decompose the circu- 
lar rotation curve measured from the high-resolution [C n] line 
into multiple components. In Sect. 5 we discuss the spatial dis- 
tribution and extent of the ISM, the surface density of the gas 
and the star formation, the ionized and molecular gas kinemat- 
ics, the gas outflow, and the dynamics of the quasar-host sys- 
tem. In Sect. 6 we summarize our results. Finally, in Appendices 
A-C, we present the channel maps of [C u] and CO (9-8) lines 
and describe the 2D Sérsic function and the tilted ring model. 
Throughout the paper we adopt a ACDM cosmology with Ho — 
67.8 km/s/Mpc, Qm = 0.3089, and €4 = 0.6911 (Planck Collab- 
oration et al. 2016). Under this cosmological assumption and at 
z = 6.0031, 1” on the sky corresponds to a physical size of 5.84 
kpc, the luminosity distance is Dg = 59.0 Gpc, and the age was 
0.9297 Gyr since the Big Bang. 


2. ALMA observations and data reduction 


We conducted ALMA band-6 observations of the [C n] line (“rest 
= 1900.5369 GHz; redshifted to voy, = 271.3851 GHz), along 
with band-4 observations of the CO (9-8) line (Vrest = 1036.9124 


GHz; redshifted to vobs = 148.0648 GHz) and the OH" (11—01) 
line (Yrest = 1033.0582 GHz; redshifted to vop, = 147.5144 GHz; 
hereafter OH+), toward J2310+1855 at z = 6.0031 from 2019 
August 08 to 20 (PI: Yali Shao; Project code: 2018.1.00597.S). 
We used 40—47 12-m antennas in the C43-7 configuration with 
a maximum projected baseline of 3.6 km for observations with 
both bands. We centered one of the 2 GHz spectral windows 
(four in total) on the redshifted [C m]/CO (9-8) line frequency, 
and used the rest of the spectral windows to observe the dust 
continuum. The total observing times are 1.60 and 1.44 hours, 
resulting in on-source integration times of 0.66 and 0.80 hours 
for the [C n] and CO (9-8) /OH* observations, respectively. We 
established the flux density scale using scans of the standard 
ALMA calibrator J2253+1608. The flux calibration uncertain- 
ties are ~5—10% for ALMA Cycle 6 bands-4 and -6 observations 
(Warmels et al. 2018). We checked the phase and water vapor by 
observing nearby calibrators of J2316+1618 and J2307+1450. 
The data were calibrated using the ALMA standard pipeline 
using the Common Astronomy Software Application (CASA'). 
The original channel width is 15.625 MHz, corresponding to 
~17 and 32 km/s for the band-6 and -4 observations, respec- 
tively, which can sample the intrinsic line widths well. The un- 
derlying dust continuum emission was subtracted in the uv-plane 
for both data sets. 


In order to improve the uv coverage and thus the final im- 
age and spectrum sensitivity, we also included archival data from 
projects — 2011.0.00206.S (PI: Ran Wang), 2015.1.00997.S (PI: 
Roberto Maiolino), and 2015.1.01265.S (PI: Ran Wang), which 
observed the [Cn] line in the 12-m array and 12-m + 7-m ar- 
rays, and the CO (9-8) /OH* (1;—0,) line in the 12-m array, to- 
ward J2310--1855, respectively. We only consider the frequency 
range overlapping with our science goals. As we used ALMA 
Cycle 0 observations for which the data weights are not cor- 
rect, before the final combination of multi-epoch data, we re- 
weighted the calibrated target data with the STATWT task in 
CASA, which attempts to assess the sensitivity per visibility and 
adjust the weights accordingly with line-free data. Finally, we 
made the line data cube from the combined calibrated data using 
the TCLEAN task in CASA with robust weighting (robust = 0.5) 
for the [C n] and CO (9-8) line in order to optimize the sensitiv- 
ity per frequency bin and the resolution of the final maps, and 
natural weighting (robust — 2.0) for the OH* line in order to im- 
prove the sensitivity. For the continuum images, we used robust 
weighting (robust — 0.5), which allows us to compare the emis- 
sion of the continuum with that of the covered lines ([C n] and 
CO (9—8)) at similar angular resolution. As the OH* line is very 
weak, we TCLEAN the data cube deeply with a threshold of 1c 
using a 1" x 1" square centered on the quasar position. For the 
rest of the lines and continuum, we TCLEAN to a level of 3c. In 
addition, during TCLEAN for the 12-m and 7-m combined data, 
we used mosaic gridder mode, which can correctly image data 
with different antenna sizes. The synthesized beam size of the 
final [C n] and CO (9-8) images are 0111 x 07092 and 07187 
x 07153, corresponding to 0.65 kpc x 0.54 kpc and 1.09 kpc 
X 0.89 kpc, respectively, at the quasar redshift. The noise lev- 
els in a 15.625 MHz channel are 0.17 and 0.10 mJy/beam for 
the [C u] and CO (9-8) lines, respectively. The root mean square 
(rms) noise in the continuum maps at 262 and 147 GHz are 0.02 
and 0.01 mJy/beam, respectively. 


! https://casa.nrao.edu/ 
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Table 1. Measurements from ALMA observations. 


[Cu] CO (9-8) OH* 
Line Intensity Map 
Weighting 0.5 0.5 2.0 
SiZebeam (4) 0.111 x 0.092 0.187 x 0.153 0.243 x 0.198 
Line Spectrum 
Z 6.0032 + 0.0001^ 6.0028 + 0.0002“ - 
FWHM (km/s) 384 + 5" 388 + 16" 225 + 117°; 240 + 744; 330 + 198° 
VS y secum (Jy km/s) 6.50 + 0.174; 6.31 +0.18? 1.44 € 0.114; 1.39 € 0.11^ 0.072 + 0.049°; 0.072 + 0.0234; —0.065 + 0.044° 
Lspectrum (10? Ls) 6.39 + 0.17%; 6.21 € 0.18^ 0.77 + 0.067; 0.74 + 0.06^ — 0.038 + 0.026°; 0.038 + 0.0134; —0.035 + 0.023* 


Continuum Map 
e 
(GHz) 


0.113 x 0.092 
262 


SiZ€beam_con 


Ycon 


0.192 x 0.156 $ 


147 - 


Notes. The parameters for the line spectra are measured using a single and double Gaussian fit for [C 1r] and CO (9-8) lines, and a triple-Gaussian 
fit for the OH+ line. The adopted [C u], CO (9-8) and OH* line spectra are shown in the right panels of Figs. 1, 2 and 3. For the [C i] and CO (9-8) 
lines: ^single Gaussian fit results; ^double-Gaussian fit results. For the OH* line: ^blue part emission component; ^red part emission component; 


*absorption component. 


3. Results 


The [C ir] and CO (9-8) lines and their underlying dust contin- 
uum are all spatially resolved. The intensity peak of OH* (11—01) 
is detected at >50 significance. We list the observational results 
in Table 1. 


3.1. The [C u] line 


The velocity-integrated intensity map, the intensity-weighted ve- 
locity map and the spectrum of the [C rr] line are presented in Fig. 
; 1. A clear velocity gradient can be seen in the velocity map. In 

addition, the [C ir] emission peak moves in a circular, counter- 
clockwise path from 219 to —-281 km/s shown in the [C u] line 
channel maps (Fig. A.1). These are the main characteristics of a 
disk with rotating gas. Similar rotating disks have been widely 
detected in the local and high-z Universe in both quasar hosts and 
galaxies without an active galactic nucleus (AGN;Wang et al. 
2013; Lucero et al. 2015; Shao et al. 2017; Jones et al. 2017, 
2021; Banerji et al. 2021; Alonso-Herrero et al. 2018; Bewketu 
Belete et al. 2021). We measured the line spectrum by integrat- 
ing the intensity of each channel from the [C n] line data cube, 
including pixels determined in the line-emitting region above 2c 
in the [C m] intensity map. The line spectrum with original spec- 
tral resolution (15.625 MHZ) is shown in the right panel of Fig. 
1, which reveals that this line has an asymmetric profile with an 
enhancement to positive velocities. 


The single Gaussian fit to the [C n] spectrum shows that the 
line width and the source redshift are consistent with our previ- 
ous Cycle 0 observations (393 + 21 km/s and 6.0031 + 0.0002, 
respectively; Wang et al. 2013). The [Cu] line flux calculated 
from the Gaussian fit to the line spectrum agrees with our previ- 
ous ALMA observations at 0’’7 resolution (8.83 + 0.44 Jy km/s; 
Wang et al. 2013) within the ~ 15% calibration uncertainty. Our 
double-Gaussian fit results are consistent with that of the single 
Gaussian fit. We also measured a consistent value by modeling 
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the observed intensity map with the 2D elliptical Sérsic model 
as described in Sect. 4.1 and shown in Table 2. 


3.2. The CO (9-8), OH" lines and the possible detection of 
para-H,O ground-state emission line 


The velocity-integrated intensity map, the intensity-weighted ve- 
locity map and the spectrum of the CO (9—8) line are presented 
in Fig. 2. A similar asymmetric component at ~—50 km/s is seen 
as in the [C n] line. 


The Gaussian fit of the CO (9-8) line reveals a line width 
and source redshift that are consistent with the previous ALMA 
Cycle 3 observations (376+ 18 km/s and 6.0031 +0.0002, respec- 
tively; Li et al. 2020b). The CO (9-8) line flux calculated from 
the Gaussian fit on the line spectrum agrees with the previous 
ALMA observations at ~ 0”7 resolution (1.31 +0.06 Jy km/s; Li 
et al. 2020b). We also obtained a consistent value by modeling 
the observed intensity map with the 2D elliptical Sérsic model 
shown in Table 2 (Sect. 4.1). 

The OH" velocity-integrated intensity map and spectrum are 
shown in Fig. 3. The peak flux density is 0.066 + 0.012 Jy/beam 
km/s (>50) in the intensity map averaged only over the emission 
component. There appears a double-horn emission profile and an 
absorption at ~—400 km/s. We used a triple-Gaussian model (red 
line in Fig. 3) to fit the spectrum, and we present the flux and line 
parameters in Table 1. The P-Cygni profile (with an absorption 
at —384 + 128 km/s from the Gaussian fit) of the OH* spectrum 
may indicate an outflow. 

Our new ALMA band-4 observations serendipitously par- 
tially cover the para-H2O (111-000) line at vrest = 1113.3430 GHz 
(redshifted to 158.9786 GHz), which is detected in emission in 
J2310+1855. The spectrum is shown in Fig. 4, where we over- 
laid a spectrum of the para-H5O (202-111) line from the ALMA 
archive. The para-H5O (111-000) line is usually detected as an ab- 
sorption line in most galaxies (e.g., Weiß et al. 2010; Yang et al. 
2013); however it shows up as a conspicuous but weak emission 
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Fig. 1. ALMA observed [C n] line. Left panel: [C n] velocity-integrated map. The white plus sign is the Hubble Space Telescope (HST) quasar 
position (RA of 23:10:38.90; Dec of 18:55:19.85). The shape of the restoring beam with a FWHM size of 0/111 x 0/092 is plotted in the bottom- 
left corner. The contour levels are [—3, 3, 6, 9, 12, 15] x rmspc;j, where rms(c,; = 0.035 Jy/beam km/s. Middle panel: [C 1] velocity map produced 
using pixels above 4c in the channel maps. The black plus sign is the HST quasar position. A clear velocity gradient can be seen. Right panel: 
[C 1] line spectrum (black histogram) overplotted with the best-fit Gaussian profiles. The spectrum is extracted from the 2c region of the source 
emitting area in the [C n] intensity map. The solid blue line presents the best-fit single Gaussian model. The solid red line represents the best-fit 
double-Gaussian model, and the dashed and dotted red lines are for each of the double Gaussians. The kinematic local standard of rest (LSRK) 
velocity scale is relative to the [C rr] redshift from our ALMA Cycle 0 observations (Wang et al. 2013). The spectral resolution is 15.625 MHz, 


corresponding to 17 km/s. The rms of the spectrum measured using the line-free spectrum is 0.65 mJy and is shown as a black bar. The [Cu] 
spectral profile is asymmetric, with enhancement on the red side. The green bar indicates the velocity range of the OH* absorption. 
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Fig. 2. Similar to Fig. 1, but for the CO (9-8) line. The FWHM size of the restoring beam is 0/187 x 07153. The contour levels are [—3, 3, 6, 
9, 12, 15, 18, 21] x rms,,, where rms, = 0.020 Jy/beam km/s. The spectral resolution is 15.625 MHz, corresponding to 32 km/s. The rms of the 
spectrum is 0.31 mJy. The CO (9—8) spectral profile is also asymmetric, with enhancement on the red side. 


feature in a few galaxies (e.g., González-Alfonso et al. 2010; 
Spinoglio et al. 2012; Liu et al. 2017). The cold (Taust ~20-30 K) 
and widespread diffuse ISM component gives rise to the ground- 
state para-H2O line (Weiß et al. 2010; Liu et al. 2017). Liu et al. 
(2017) find the ratio of para-H5O (111-000) and para-H20 (2 9- 
111) to be < 1 in galaxies in which both lines are detected as 
emission lines. As the ground-state para-H2O line is very close 
to the edge of the spectral window in our observations, we only 
detect it partially and may not cover the peak of its spectrum. 
Thus, we only simply report the detection of this line and will 
not discuss it further in what follows. 


3.3. The dust continuum emission 


The dust continuum distribution, overplotted with that of the 
[C u] and CO (9-8) lines, is shown in the left and middle panels 
of Fig. 5. We study the brightness distribution and the effective 


radius (half-light radius) for both lines and their underlying dust 
continuum emission in Sect. 4.1. The measured intensities are 
shown in Table 2. The z ~ 6 source is spatially resolved with our 
data, so we investigate the brightness distributions with 2D ellip- 
tical Gaussian functions. However, for the low-resolution data, 
the observed brightness distribution is dominated by the large 
beam, so the CASA IMFIT tool, which fits a 2D elliptical Gaus- 
sian, can be used effectively. The flux densities (listed in Table 
2) of the [C n] and CO (9-8) underlying dust continua measured 
using the fits in Sect. 4.1 are consistent with the ALMA Cycle 
0 data (8.91 + 0.08 mJy; Wang et al. 2013) and ALMA Cycle 3 
data (1.59 + 0.04 mJy; Li et al. 2020b), respectively, measured 
using the CASA IMFIT tool including additional calibration un- 
certainties. As shown in Table 2, the effective radii of the [C n] 
and CO (9-8) lines are larger than that of their underlying dust 
continuum emission, as we discuss in Sect. 5.1.1. 
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spectrum (black histogram) overplotted with the best-fit triple-Gaussian profile (red line). The dashed and dotted blue lines present the blueward 
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Redshift 


5.98 6.00 6.02 
Ow UU TXELPLTÀ 


| P-H20(11i — 000) 
| p-H50(2o; — 111) x0.28 


0.3— 


0.4 


0.2 
0.1 


0.0 ——— FE 


rit 


-1500 -1000 -500 (0 500 
LSRK velocity (km/s) 


el eel ei ey ae Ds i 
1000 1500 


Fig. 4. Spectra of the para-H2O (111-000) line (black histogram) and 
the scaled para-H5O (202-111) line (blue histogram). The LSRK velocity 
scale is relative to the [C rr] redshift from our ALMA Cycle 0 observa- 
tions (Wang et al. 2013). The spectral resolution for both lines is 15.625 
MHz, corresponding to ~32 km/s. The observations may only cover the 
red tail of the para-H5O (111-000) emission line. 
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4. Analysis 
4.1. Image decomposition 


In order to investigate the ISM distribution, we used a 2D ellip- 
tical Sérsic function (Eq. B.3) to reproduce the observed inten- 
sity maps of the [Cu] and CO (9-8) lines and their underlying 
dust continuum. The 2D elliptical models represent the inclined 
brightness distribution. We first convolved the model with the 
restoring beam kernel, and then determined the best fit Sérsic 
model by comparing the convolved model and the data pixel-by- 
pixel using Markov chain Monte Carlo (MCMC) method with 
the emcee” package (Foreman-Mackey et al. 2013). The best-fit 
models are shown in the left panels of Figs. 6 (for the [C u] line), 
7 (for the CO (9-8) line), B.2 (for the dust continuum underly- 
ing the [Cn] emission), and B.3 (for the dust continuum under- 
lying the CO (9-8) emission). During the fitting, we considered 
uniform error maps as we only used the central regions of the 
Observed intensity maps. Motivated by the possible existence of 
a torus surrounding the AGN in the AGN unification model and 
from observations of several molecular lines and dust contin- 
uum (e.g., Hónig 2019; Combes et al. 2019; García-Burillo et al. 
2021), we also add an additional unresolved nuclear component 
(which is concentric with the extended component) to represent 
the dusty and molecular torus for the CO (9—8) line and the dust 
continuum emission. For the [C n] line, we do not include this 
component, as the observed [C ri] emission is suppressed by the 
high dust opacity, as we discuss in Sect. 5.1.1. We note that in 
the following analysis and discussion on the spatially resolved 
gas and dust content in the quasar host galaxy, we remove the 
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Fig. 5. Dust continuum images at different frequencies and their radial profile comparison. Left and right panels: [C ir] and CO (9-8) underlying 
dust continuum maps at 262 and 147 GHz, respectively. The black plus sign is the quasar position from HST snapshot observations. The bottom- 
left ellipse in each panel shows the restoring beam with FWHM sizes of 07113 x 07092 and 0/192 x 0/156, respectively. The black contours 
are [—3, 3, 9, 27, 81] x rms, for both maps. The rms,,., at 262 and 147 GHz are 0.02 and 0.01 mJy/beam, respectively. The overplotted white 
contours are from [C n] and CO (9-8) emission lines, and the contour levels are the same as those of the black contours in the left panels of Figs. 1 
and 2. Right panel: Comparison of the surface brightness of the continuum at two different frequencies. The blue and red symbols with error bars 
are for the scaled 262 and 147 GHz dust continuum, respectively. The error bar represents the deviation of the values of all pixels in each ring used 


to do the aperture photometry. 


emission associated with the dusty and molecular AGN torus 
(i.e., subtract the unresolved nuclear component model from the 
observed image). To determine which model (Sérsic or point 
source+Sérsic) can better fit the line and/or the dust continuum 
* data, we first measured the pixel-to-pixel rms within a 1” x 1” 
square at the center of the residual map (the measured values are 
shown in the captions of Figs. 7, B.2 and B.3), and then calculate 
the distribution of residual values inside the square (the top-right 
histogram of Figs. 7, B.2, and B.3). In addition, we compared 
the surface brightness profiles of both the modeled and observed 
intensity maps with projected rings (corrected by the inclination 
angle and position angle given in Table 3 determined from the 
kinematic models in Sect. 4.3). 


The parameters of the best-fit image decomposition for each 
emission line and its associated continuum are listed in Table 2. 
For the [Cn] line (shown in Fig. 6), the best-fit Sérsic index is 
0.59, which is smaller than a typical exponential disk with a Sér- 
^ sic index of 1. It is also smaller than that of the CO (9-8) line and 
the dust continuum of our quasar J2310+1855. This discrepancy 
may be due to the high dust opacity, rather than the intrinsic 
distribution of the [C i] line, as we discuss in Sect. 5.1.1. For 
the CO (9-8) line, the best-fit Sérsic index is 2.01 in the case of 
a single Sérsic model. When we add a point component in the 
center, the Sérsic index of the extended component decreases to 
1.21. The statistics on the residual maps for the one- and two- 
component scenarios appear identical. But the two-component 
model better matches the observed CO (9-8) within the central 
^] kpc region, as shown by the surface brightness difference dis- 
tribution in Fig. 7. The [C ii] and CO (9-8) underlying dust con- 
tinua, similar to the CO (9-8) line, seem better described by a 
combination of a point component and an extended Sérsic com- 
ponent, as shown in Figs. B.2 and B.3. And the Sérsic index 
for the dust continuum emission is a little bit larger than 1. The 
dust continuum is more concentrated than the gas emission, as 
measured by the half-light radius. This is consistent with what 
is found in observational and theoretical studies of high-redshift 
galaxies (e.g., Strandet et al. 2017; Tadaki et al. 2018; Cochrane 
et al. 2019). We discuss the spatial distribution and extent of the 


ISM in Sect. 5.1. The central positions (see Columns 3-4 in Ta- 
ble 2) of the gas and the dust continuum are identical. 


4.2. Dust diagnostic 


With the dust emission size from the image decomposition, and 
following Weiß et al. (2007) and Walter et al. (2022), we are 
able to constrain the dust temperature and the dust mass with 
a general gray-body formula instead of using the optically thin 
approximation as done in our previous work (Shao et al. 2019): 


S, E Opp [By(Taust) = By,(T cwp)ll1 e» exp(-7,)](1 * a, (1) 
and 
T= Ko] Viet Y Maust, app / (D4 Qapp)s (2) 


where S, is the observed flux density for a target at redshift z. 
By(Taust.) and B,(Tcwp) are the black-body functions with dust 
temperature Tgys and the cosmic microwave background (CMB) 
temperature Tcmp, respectively. Qapp is the apparent solid angle. 
In principle, Qapp can be different for different wavelengths. Due 
to the dust continuum sizes remaining constant within ~20% at 
observed frame wavelengths from 500 um to 2 mm for z ~ 1 — 5 
main-sequence galaxies from simulations (Popping et al. 2022), 
and our findings for the similar dust sizes (within ~10%) at 
wavelengths of 158 and 290 um in the rest frame, we do not 
consider changes in this quantity below. T, is the optical depth, 
which is a function of frequency v (rest-frame) and dust mass 
surface density Maust, app/ (DZ Qapp). Maust, app is the apparent dust 
mass. Da is the angular distance. 6 is the emissivity index. We 
here adopted the absorption coefficient per unit dust mass «o of 
13.9 cm?/g at the reference frequency vref of 2141 GHz (Draine 
2003; Walter et al. 2022). 

Motivated by the AGN contribution to the strength of the 
dust continuum emission of the quasar-host system (e.g., Di 
Mascia et al. 2021), we decomposed the UV-to-FIR SED of 
J2310+1855 shown in Fig. 8. With the best-fitting UV/optical 
power law (red lines) and NIR to mid-infrared (MIR) CAT3D 
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Fig. 6. Image decomposition of the [C1] line. Left and middle panels: 2D elliptical Sérsic model and the residual between the observed and 
modeled [C i] intensity maps. The shape of the [Cm] synthesized beam with a FWHM size of 0111 x 07092 is plotted in the bottom-left corner 
of each panel. We measured the residual rms to be 0.037 Jy/beam km/s inside the dashed gray square with a side length of 1". Top-right panel: 
[Cn] luminosity surface brightness (black squares and red diamonds with error bars are measurements from the observed intensity map and the 
2D elliptical Sérsic model, respectively) at different radii, measured using elliptical rings with the ring width along the major axis half (005) that 


< , ofthe [C n] clean beam size, the rotation angle equal to PA (=199°), and the ratio of semiminor and semimajor axis — b/a of cos(i) (i = 42°), where 


PA and i come from the [C n] line kinematic modeling. The error bar represents the standard deviation of the values of all pixels in each ring. The 
solid and dashed gray lines are three times the [Cm] luminosity surface brightness limit measured in the emission-free region with an elliptical 
annulus that is the same with the ones used to measure the [C n] luminosity surface brightness but with a larger radius, and its corresponding rms. 


=. Bottom-right panel: Luminosity surface brightness difference (red diamonds) between the measurements from the observed intensity map and the 


2D elliptical Sérsic model. Note that the vertical scale is linear and in units of 107 Lo /kpc?. 


AGN torus model (brown lines; Hónig & Kishimoto 2017) in 
Shao et al. (2019), and the source size of twice the effective ra- 
dius of the dust continuum underlying CO (9-8) listed in Table 2, 
we derived an average dust temperature Taust of 5233 K, a total 


dust mass Maust, app Of 127200 x 10? M. and an overall emis- 


sivity index B of 1.00055. Our dust temperature is smaller than 
the value of 71 + 4 K derived by Tripodi et al. (2022), and our 


‘s dust mass is ~3 times higher. This is due to the differences in 


the AGN dust torus model used, which has a significant contri- 
= bution in the rest frame wavelength range <100 um in our case, 
whereas the impact of the AGN dust torus on the total dust con- 
tinuum emission in Tripodi et al. (2022) is limited x50 um. As 
a result, the peak of our gray-body model is at a longer wave- 
length, which leads to a lower dust temperature and a higher dust 
mass according to Eqs. 1 and 2, as well as a smaller gas-to-dust 
ratio (GDR) compared with the dust property results in Tripodi 


et al. (2022). The FIR and IR luminosities are 8.847245 x 10? Lo 


and 11097 x10! Lo, respectively, by integrating the gray-body 
model (which can be taken as purely star-forming heated dust 
emission) from 42.5—122.5 um and from 8—1000 um. The aver- 
age SER is 2055% M,/yr calculated from the above-mentioned 
IR luminosity and the formula in Kennicutt (1998). 

With the resolved dust continuum at two different frequen- 
cies, we are able to investigate the radial distributions of the dust 
temperature, dust mass and SFR. We first used circular annuli to 
do aperture photometry on both the [C i1] and CO (9-8) underly- 
ing dust continuum maps (after removing the central point com- 
ponents that may come from the AGN dust torus) with ring width 
of 0” 1. The right panel of Fig. 5 shows the comparison of the sur- 
face brightness of the continuum at observed frame 262 and 147 
GHz. The slight differences in the dust continuum profiles at dif- 


ferent wavelengths may indicate different distribution of the dust 
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temperature along with distance. For each ring we fitted the mea- 
sured flux densities with Eqs. 1 and 2 using the £ value of 1.90 
from our SED fitting, to get the average dust temperature, dust 
mass and dust optical depth at each radius. Finally, with fixed 
gray-body formula for each ring, we calculated the IR luminos- 
ity and SFR. The results are shown in Fig. 9. The dust tempera- 
ture drops with increasing radius as seen in the top-left panel. A 
similar dependence is found in a star-forming galaxy at z — 7.13 
(Akins et al. 2022). We fit a power law, Taust « p 0415020 The 
radial dependence of X4ust mass and Xsrg are consistent with an 
exponential distribution X(r) = Xo exp”! ^, where r, is the ex- 
ponential scale length and Xo is the central surface density. The 
best-fit relation for Zdust_mass is shown in the top-right panel of 
Fig. 9 with r,aus; mass = 0.77 + 0.27 kpc, corresponding to an 
effective radius of 1.29 + 0.45 kpc. It is about two times larger 
than the dust continuum effective radius (i.e., ~0.6 kpc) shown in 
Table 2. This may indicate that the single-band dust continuum 
emission is not a good tracer of the dust mass, which is consis- 
tent with the finding from numerical simulations (Popping et al. 
2022), and that the dust emission depends more strongly on dust 
temperature than on dust mass. Thus, one requires at least two 
bands of high-resolution imaging to map the dust temperature 
across the galaxy disks when using the dust continuum emission 
of galaxies as a reliable tracer of the dust mass distribution. The 
total dust mass summing over all rings is (1.84 + 0.69) x 10° Mo, 
which is consistent with the value measured from the SED fit- 
ting. An exponential fit for Xspr is shown in the bottom-right 
panel of Fig. 9 with rs sgg = 0.38 + 0.09 kpc, corresponding to 
an effective radius of 0.64+0.15 kpc. It is consistent with the dust 
effective radius, which means that the resolved dust continuum 
emission is very closely linked to the SFR distribution. These 
are consistent with the conclusions from TNGSO star-forming 
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; Fig.7. Image decomposition of the CO (9-8) line. Top-left panels: 2D elliptical Sérsic model and the residual between the observed and modeled 


CO (9-8) intensity maps, from left to right. We measured the residual rms inside the dashed gray square with a side length of 1". Bottom-left 
panels: Same as in the top panels but for the Point+2D elliptical Sérsic modeling. The peak values and rms within the dashed gray square are 0.061 
and 0.021, and 0.059 and 0.020 Jy/beam km/s for the Sérsic and P+Sérsic residual maps, respectively. The shape of the CO (9-8) synthesized 
beam with a FWHM size of 07187 x 07153 is plotted in the bottom-left corner of each panel. Top-right panel: Distributions of the residuals 
for pixels inside the dashed gray squares (Sérsic modeling: blue histogram; P+Sérsic modeling: red histogram). Middle-right panel: CO (9-8) 
luminosity surface brightness (black squares, blue circles, and red diamonds with error bars are measurements from the observed intensity map, 
the 2D elliptical Sérsic model, and the Point+2D elliptical Sérsic model, respectively) at different radii, measured using elliptical rings with the 
ring width along the major axis half (071) that of the CO (9-8) clean beam size, the rotation angle equal to PA (=198°), and the ratio of semiminor 


and semimajor axis — b/a of cos(i) G = 43°), where PA and i come from the CO (9—8) line kinematic modeling (listed in Table 3). The error bar 
represents the deviation of the values of all pixels in each ring. The solid and dashed gray lines are three times the CO (9-8) luminosity surface 


! brightness limit measured in the emission-free region with an elliptical annulus that is the same with the ones used to measure the CO (9-8) 


luminosity surface brightness but has a larger radius, and its corresponding rms. Bottom-right panel: Luminosity surface brightness difference 
between the measurements from the observed intensity map and the ones from the 2D elliptical Sérsic model (blue circles) and from the Point+2D 
elliptical Sérsic model (red diamonds). Note that the vertical scale is linear and in units of 10°L,/kpc?. 


4.3. Gas kinematic modeling 


The [Cu] and CO (9-8) lines can be used to trace the kine- 
matic properties of atomic, ionized, and molecular gas in the 
quasar host galaxy. As shown in Fig. A.1, the [Ci] emission 
peak moves in a circular path with frequency. As shown in Figs. 
10 and C.2, the velocity fields traced by both [C ir] and CO (9-8) 
show clear velocity gradients. In addition, the position-velocity 
diagrams along the kinematic major axes have an "S" shape (es- 
pecially apparent for the [Cu] line; however, the resolution of 
the CO data is roughly two times lower than that of the [Cu] 
data, so the "S" shape is not obvious). These are consistent with 
a rotating gas disk. We are able to construct a 3D disk model for 
these data cubes with the package called 3D-Based Analysis of 
Rotating Object via Line Observations ??BaroLo*; Di Teodoro 
& Fraternali 2015). 

The ?PBanoro software fits 3D tilted ring models to spec- 
troscopic data cubes. For each ring, the algorithm builds a 3D 


galaxy simulations that the single-band dust emission is a less 
robust tracer of the dust distribution, but is a decent tracer of 
the obscured star formation activity in galaxies (Popping et al. 
2022). The total SFR summing over all rings is 341545746 Mo/yt, 
which is consistent with the value measured from the SED fit- 
ting. We suggest adopting the total SFR derived from the SED 
fitting, as the SFR value for each ring is measured from only two 
data points at different wavelengths and thus has a lot of uncer- 
tainty. As shown in the bottom-left panel of Fig. 9, as the dust 
mass surface density decreases with the galactocentric radius, 
the dust emission in both bands becomes less optically thick. 
The rest-frame wavelength of both dust bands (~290 um for the 
CO (9-8) underlying dust continuum, and ~158 um for the [C 1n] 
underlying dust continuum) is shortward of the Rayleigh-Jeans 
tail (i.e., a rest-frame wavelength of ~350 um). The dust mass 
surface density is very high (i.e., Xaus mass > 108 M; /kpc?), and 
the optical depth for the dust continuum emission at both fre- 
quencies is above 0.1 inside ~1 kpc. 


3 http: //editeodoro. github.io/Bbarolo/ 
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Fig. 8. SED decomposition toward J2310--1855. Left panel: Rest-frame UV-to-FIR SED fitting. The red points with error bars or downward arrows 
are observed data (see Table 3 in Shao et al. 2019 for details). The pink lines represent the UV/optical power law from the accretion disk. The 
brown lines correspond to the CAT3D AGN torus model (Hónig & Kishimoto 2017). The physical properties derived for these two components 
are detailed in Table 4 in Shao et al. (2019). The green lines correspond to a gray-body profile (Eqs. 1 and 2) associated with star formation 
activity in the quasar host galaxy. The black lines are the sum of all components. The fit employed the MCMC method with the emcee package 
(Foreman-Mackey et al. 2013). We visualized the model uncertainties with shaded areas by randomly selecting 100 models from the parameter 
| space. Right panel: Corner map for the parameters (dust temperature, dust mass, and dust emissivity index) associated with the gray-body model 

for the FIR region SED fitting (green lines in the left panel). The contours are drawn at 1 — exp(-n? /2) (m = 0.5,1,1.5,2) of the volume. The 


vertical dashed lines show 16th, 50th (median), and 84th percentiles. 


disk model of the gas distribution in both the spatial and veloc- 
| ity axes, and then convolves it with the restoring beam of the 
observed data cube. Finally, it compares the convolved data set 
with the observed one. The geometry of the tilted ring model 
can be seen in Fig. C.1 and the corresponding geometrical and 
kinematic parameters (e.g., the inclination angle i; the position 
angle ¢; the rotation velocity V;4,) are described in Appendix C. 
; The 3D tilted ring model performs better than the standard 2D 
modeling on the velocity fields. For example, a common prob- 
lem when deriving the kinematic properties from the velocity 
fields is the beam smearing effect, which smears the steep ve- 
locity gradient especially in the central region. And there exists 
differences among the velocity fields measured using different 
methods (e.g., an intensity-weighted velocity field from CASA, 
versus a mean gas velocity map based on Gaussian fitting us- 
ing AIPS^. The ?PBanoro algorithm avoids the beam smearing 
with the convolution step and directly conducts the modeling 
on the data cube. In addition to the rotation velocity for each 
ring, "PBARoLO can also measure the asymmetric-drift correc- 
tion (VA), which is caused by the random motions and can be 
directly measured with the velocity dispersion and the density 
profile (e.g., Iorio et al. 2017). There is no direct way to estimate 
parameter errors, and "PBARoLOo adopts a Monte Carlo method: 
it calculates models by changing the parameter values with ran- 
dom Gaussian draws centered on the minimum of the function, 
once the minimization algorithm has converged. 

We took the ring width W,ing as half of the clean beam size 
(i.e., Nyquist sampling): 005 and 0"1 for the modeling of the 
[C u] and CO (9-8) lines, respectively. The central radius of each 
ring is Wis, X t where ¢ is the t-th ring. Eight rings are adopted 


^ http://www.aips.nrao.edu/ 
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for the kinematic modeling on our high-resolution [Cn] data. 
The [Cn] emission extends to a radius of ~2.8 kpc as shown by 
the [Cn] surface brightness distribution (top-right panel in Fig. 
6). We chose to use data within a radius of ~2.4 kpc, taking ad- 
vantage of higher S/N data. We should note that as the projected 
[C u] kinematic minor axis (~2.5 kpc) is longer than the kine- 
matic major axis, we will lack some information near the kine- 
matic major axis for the outer two rings. However, since we have 
plenty of information from near the kinematic minor axis and in 
the region between these two axes, we still have enough data to 
carry out the kinematic modeling in the outer two rings. The ini- 
tial guess of the position angle and the inclination angle are 200° 
and 40°. The initial position angle are roughly measured from 
the velocity maps based on the definition of the position angle 
(see Appendix C). One method for roughly determining the in- 
clination angle, assuming an intrinsic round disk, uses the photo- 
metric minor and major source size ratio (size): i = cos" (Fsize). 
As the kinematic major axis of the [C 11] line is far from its pho- 
tometric major axis and [Cu] kinematic major axis is shorter 
than its kinematic minor axis, but the kinematic major axis of 
the CO (9-8) line is close to its photometric major axis (see Ta- 
ble 2) and CO (9-8) kinematic major axis is longer than its kine- 
matic minor axis, we derived the initial inclination angle from 
the CO (9-8) minor and major source size ratio. In addition, dur- 
ing the fitting, we adopt pixel-by-pixel normalization, which al- 
lows the code to exclude the parameter of the surface density 
of the gas from the fit. It means that we force the value of each 
spatial pixel along the spectral dimension in the model to equal 
that in the observations, which allows a non-axisymmetric model 
in density and avoids untypical regions (Lelli et al. 2012). We 
present the fitted parameters in Table 3. The average inclination 
angle and position angle derived from the [Cn] and CO (9-8) 
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; Fig. 9. Radial distributions of the dust temperature (top left), dust mass surface density (top right), optical depth (bottom left) at [C n] and CO (9-8) 
underlying dust continuum frequencies (red dots and blue squares with error bars, respectively), and the SFR surface density (bottom right). The 
purple, gold, and green lines are best-fit lines. Note that the uncertainties of the optical depth at different wavelengths are same in the bottom-left 
panel but appear different due to different ranges of the left and right axes. The error bars of the SFR surface density (which is associated with 
the dust temperature and dust mass) in the bottom-right panel are measured with a Monte Carlo method: we calculate samples by changing the 
parameter values with random Gaussian draws centered on their best-fit values and deviated by their errors shown in the top-left and top-right 
panels, and the lower and upper values are taken as the 16th and 84th percentiles. 


lines are consistent. We plot the intensity, velocity, velocity dis- 
persion, and position-velocity maps and the line spectra mea- 
sured from the modeled line data cubes in Figs. 10 and C.2 for 
the [C rr] and CO (9-8) lines, respectively. The rotation velocity, 
gas velocity dispersion, and asymmetric-drift correction for each 
ring are shown in Fig. 11. The circular rotation velocities (the 
rotation velocities corrected by the asymmetric drifts; see Eq. 3) 
presented in Fig. 12 for [C n] and CO (9-8) lines are consistent 
with each other. 


4.4. Rotation curve decomposition 


In order to quantify the dynamical contribution of each matter 
component, we performed a decomposition of the circular rota- 
tion curve measured from the high-resolution [C n] line. 


The circular velocity (V.) directly traces the galactic grav- 
itational potential () and can be measured by correcting the 
rotation curve for the asymmetric drift: 


O® 
R OR (3) 
where Vis is the rotation velocity, Va is the asymmetric drift 
correction caused by random motions and can be modeled given 
the velocity dispersion and the density profile (e.g., Iorio et al. 
2017) by ?PBanoLo, and 6 is the sum of the potentials of the 
different mass components. 

We consider four matter components (black hole, stellar, gas, 
and dark matter) that contribute to the total gravitational poten- 
tial of the quasar-host system. Thus, the circular velocity can be 
expressed as 


-du 2 2 2 
Vo = WV eH + Vitor + Vias + Vp: 


2 2 2 
= Ve E Viot * Vi 


(4) 
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Fig. 10. Observed and modeled [C n] line of J2310+1855. Upper-left panels: Intensity, velocity, and velocity dispersion maps for the observed 
data from top to bottom, which are labeled “D.” Upper-central panels: Same as the left panels but for the best-fit modeling from ?PBAROLO, 
labeled “M.” Upper-right panels: Same as the left panels but for the residuals between the observed data (“D”) and the modeled ones (“M”), which 
are labeled “R” (note that the dynamical range of these panels are different from others). The black cross in each panel marks the center of the 
rotating gas disk. The shape of the synthesized beam with a FWHM size of 0111 x 07092 is plotted in the bottom-left corner of each panel. In 
the velocity field panels, the dashed black lines are the kinematic major axis of the gas disk, and the plotted solid gray dots represent the ring 
positions. Lower-left and lower-middle panels: Position-velocity maps extracted along the kinematic major and minor axes from the observed data 
cube (blue contours) and the PBAROLO model cube (red contours). The contour levels are [2, 2, 6, 10] x 0.15 mJ y/beam for both the data and the 
model. Lower-right panel: [C ir] spectra extracted from the observed data cube (blue histogram) and the PBAROLO model cube (red histogram). 
The residual spectrum between the data and the model is shown as the black histogram. The spectral resolution is 62.5 MHz, corresponding to 
68 km/s. In the "PBAROLO modeling, we used a re-binned (four original channels) data cube in order to optimize the data S/N per frequency and 
velocity bins and the sampling of the FWHM of the [C rr] line. 


where Vgu, Vstar, Veas and Vpy are the contributions of the black First, the Keplerian velocity due to the central black hole — 
hole, stellar, gas and dark matter components to the circular ve- Vgy is 
locity. The following applies for each component at a radius R: 


GMpu 
R > 


(5) 


Van = 
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Fig. 11. Measured rotation velocity (left panel), gas velocity dispersion (middle panel), and asymmetric-drift velocity (right panel) as a function of 
radius (the central radius of each ring) measured from the 3PBAROLO modeling for the [C rr] (red points with error bars) and CO (9-8) (blue points 


with error bars) lines, respectively. 


where G and Mpy are the gravitational constant and the black 
hole mass, respectively. We measured a circular velocity of 
249+, km/s at the innermost ring with a central radius of 0.29 
kpc. Only considering the gravitational potential caused by the 
central black hole at the radius of 0.29 kpc, we calculated a black 
4 hole mass upper limit of 4.18*° pee x 10? Me. The black hole mass 
for our target J2310+1855 measured from the Mgn and Civ 
lines is (3.15-5.19) x10? Ms (Jiang et al. 2016), which corre- 
. sponds to a Vgg range of 216-277 km/s at R = 0.29 kpc. The 
black hole contribution to the total potential is significant in the 
innermost region, but some values of Vgy exceed the measured 
value at 0.29 kpc. We set the black hole mass as a free parameter 


^j in the dynamical modeling. The black hole sphere of influence 


| (ry) is defined as the region within which the gravitational po- 
tential of the black hole dominates over that of the surrounding 
stars, which can be measured with rp = GMgy/ 02, where Ostar 
is the velocity dispersion of the surrounding stellar population. 
When Mpy = 2.97 x 10? Mo (from this dynamical modeling) 
and Ostar = 250 km/s (see Sect. 5.5.2), ry = 0.21 kpc. When 
: Mey = 10? Mo (the order value of the black hole mass measured 
from NIR spectral lines with the local scaling relations; Jiang 
et al. 2016; Feruglio et al. 2018) and Ostar = 100 km/s (assuming 
it is equal to the maximum velocity dispersion measured from 
) our kinematic modeling of the gas), rn = 0.43 kpc. Our high- 
resolution ALMA [C ri] observations can zoom into the sphere 
of influence considering that our innermost point is at 0.29 kpc 
and, thus, can be used to derive the dynamical mass of the cen- 
tral black hole. Beyond the sphere of influence, the gravitational 
dominance of the black hole quickly vanishes, which is shown as 
the yellow line of the left and middle panels of Fig. 12. We are 
not spatially resolving the sphere of influence considering the 
ratio 75 /r;e, (i.e., 0.4—0.8 in this work) between the radius of the 
black hole sphere of influence and the spatial resolution of the 
data. However, gas dynamical studies (e.g., Ferrarese & Merritt 
2000; Graham et al. 2001; Marconi & Hunt 2003; Valluri et al. 
2004) have addressed that resolving the sphere of influence is 
an important but not sufficient factor to dynamically estimate the 
black hole mass, and the ratio rj /r;«, can be taken as a rough in- 
dicator of the quality of the black hole mass estimate. Ferrarese 
& Ford (2005) have summarized a list of black hole mass detec- 
tion based on resolved dynamical studies in their Table II, which 
shows a ry/r;e, range of 0.4—1700. Our rj/r;«, ratio is just within 
the above-mentioned range. In summary, we are able to derive a 
dynamical mass of the black hole from our [C u] data. 


Second, the stellar component can be described by a Sérsic 
profile (e.g., Terzić & Graham 2005; Rizzo et al. 2021), giving 
rise to 


Woks G Mstar VNstar(3 EM p), b(R/Re; star)! I9) 
id R I (naa (G d p)» , 


where Mstar, Re: star and Mstar are the total stellar mass, the effective 
radius and the Sérsic index. y and I are the incomplete and com- 
plete gamma functions. p and b are related to the Sérsic index by 
pe=l. = 0.6097 /Nstar + 0.05563/n2,,, (when 0.6 < n,4, < 10, 
and 10? < R/Re: star < 103) and b = 2nNgtar — 1/3 + 0.009876 /Ngtar 
(when 0.5 < star < 10). Due to the limited resolution (~600 pc) 
of our ALMA data (giving us few data points), the lack of re- 
solved rest-frame optical or NIR data to constrain the two stellar 
components (i.e., bulge and disk), and the strong degeneracies 
between the two subcomponents, following the fitting method 
from Rizzo et al. (2021), we adopted a single Sérsic element to 
globally describe the stellar component. 

Third, the gas component can be represented by a thin expo- 
nential disk (Binney & Tremaine 1987). Thus, 


(6) 


2G Maas 
Veas = VP Uo) Koy) - hpKi)]. (7) 


gas 


where M,,, and Rg4 are the total gas mass and the gas scale 
radius, /;(y) and K;;(y) are the first and second kind modified 
Bessel functions of zeroth (ii = 0) and first (ii = 1) orders, and 
y is given by y = R/(2Rg,s). The CO (9-8) image decomposi- 
tion in Sect. 4.1 suggests an approximately exponential distribu- 
tion of the molecular gas. We fixed the Sérsic index to be 1 for 
the CO (9-8) gas distribution, and got a scale radius of ~0.78 
kpc, which can be taken as the overall molecular gas scale ra- 
dius under the assumption of an exponential gas disk. With the 
detected CO (2-1) line toward J2310+1855 (vS, = 0.18+0.02 Jy 
km/s; Shao et al. 2019), and assuming Looo- n% Leoa -0) (Car- 
illi & Walter 2013), we are able to measure the gas mass with 
an assumed CO-to-gas mass conversion factor (aco), namely: 
Mgs = æcoLcoa-o The only free parameter for the gas compo- 
nent 1S Qco. 

Finally, the dark matter component contributes much less to 
the total gravitational potential than the baryonic components in 
the innermost regions. However, it is required to fit a flat rotation 
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Table 2. Parameters derived from best-fit of the image decomposition. 
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Notes. Column 1: the name of the ISM tracer. Column 2: the best-fit model. “S" and “P" present the Sérsic and point component, respectively. Columns 3-4: source position — RA and Dec, 
respectively. Column 5: the Sérsic index of the ISM (the gas and the dust emission). Columns 6-7: the scale lengths along the major and minor axes, respectively. Column 8: the rotation angle 


defined in Appendix B. Columns 9-10: the half-light radii. The values in boldface are measurements near the kinematic major axes. Columns 11-12: the line flux for the point component and 


extended component, respectively. Columns 13-14: the dust continuum flux density for the point component and extended component, respectively. Column 15: the fraction of the nuclear component 
in the whole emission. Note that the errors of these parameters are all fitting-type errors. The complex correlated noise is not included. As claimed by the task reference of the CASA IMFIT, the 


correlated noise can bring 446 effect when 5c detection of the emission for a 2D elliptical Gaussian fitting. 


curve in the outer parts in the Galaxy and other local galaxies 
(e.g., Carignan et al. 2006; Sofue et al. 2009). We consider the 
dark matter halo as a Navarro-Frenk-White (NFW; Navarro et al. 
1996) spherical halo. These assumptions lead to 


GMpy | In(1 + cx) — cx/(1 + cx)) 
Rpm x ln(1-c)-c/(1-c) ' 


VDM = (8) 
where Mpp and Rp» are the virial mass and radius (Moo) and 
R»oo in Navarro et al. 1996), respectively, c is the concentra- 
tion parameter, and x = R/R2o is the radius in units of the 
virial radius. The virial mass and radius are correlated with 


the critical density (pci): M200 = 200pc (47/3) R3, where 


Perit = 3H2/(81G) and H; = Ho YOm(1 + 23 + Qa (where z is 
the redshift). We calculate the concentration parameter (c) from 
the mass-concentration relation estimated in N-body cosmolog- 
ical simulations by Ishiyama et al. (2021). The only free param- 
eter for the dark matter component is Mpm. 

In summary, for the dynamical modeling, we have six free 
parameters (Mgu, Mstar, Mstar, Re: star, (co, Mpm) and eight data 
points (the circular velocities corrected by the asymmetric drift 
shown in Fig. 12, which have already been corrected for the in- 
clination angle from ?PBAnoro). During the fitting with the emcee 
package (Foreman-Mackey et al. 2013), a loose prior constraint 
of [10°, 10!°] Mg is adopted for Mgy. The quantities of Mstar, “star 
and Re: star are tightly coupled together, and we use prior ranges 
of [107, 10!']Mo, [0.5, 10], and [0, 2.5] kpc, respectively. As for 
Qco, we used a range of [0.2, 14] Mo/(K km/s pc’), which is 
measured from a sample of nearby AGN, ultra-luminous infrared 
galaxies, and starburst galaxies (Mashian et al. 2015). For Mpm, 
we consider a uniform prior of 10!° to 10? M, dark matter halo. 

We present the best-fit decomposition for the circular ro- 
tation curve in the left panel of Fig. 12, and the physical pa- 
rameters measured from the dynamical modeling in column 


(a) of Table 4. The black hole mass of 2.97*05! x 10? M; is 


consistent with the measurements — (4.17 + 1.02) x 10? and 
(3.92 + 0.48) x 10? M, from the Mg u and C rv lines (Jiang et al. 
2016), and (1.8 + 0.5) x 10? M, from the Civ line (Feruglio 
et al. 2018). The [C ri] emission can be taken as a molecular gas 
mass tracer in galaxies (e.g., Zanella et al. 2018; Madden et al. 
2020; Vizgan et al. 2022) with a conversion factor (ajc) by 
Meas = Q(cu]Licuj. The derived gas mass corresponds to a aqcj 
of 325 Mo/Lo. Given that there are only two degrees of free- 
dom, and the strong degeneracies among some parameters (i.e., 
Mestar, Mstar and Re: star), we Cannot constrain most of these param- 
eters well with the ~600 pc resolution data we have. We tried 
another experiment by fixing the black hole mass to a smaller 
value of 1.8 x 10? Me (Feruglio et al. 2018), and present the dy- 
namical results in the middle panel of Fig. 12 and column (b) 
of Table 4. In both situations, the derived stellar mass is on the 
order of 10? Mo (however with large uncertainty). And when fix- 
ing the black hole component, the stellar component (green line 
in the middle panel of Fig. 12) has a bump in the inner region. 
This may indicate that a massive stellar bulge already formed at 
z- 6. 


5. Discussion 
5.1. The spatial distribution and extent of the ISM 


5.1.1. Comparisons among different ISM tracers 


The different distribution and extent of various ISM tracers may 
be due to dissociation regions heated by starbursts and/or AGN 
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Table 3. Kinematic parameters derived from the 7?>BAROLO modeling. 


Data cube i PA i PA Viat Vmax O ext O med Vfat/ Oext Vinax / O med 
C) C) C) C) (km/s) (km/s) (km/s) (km/s) 

(1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) 

[Cn] Q7-45)0 9 (191-21 S 4243 199+8 25072 28343 37: 54$ TE 3 

CO(9-8) (41-43)'G (195-200) 79 4351 198-23  - 25H} - 6478 7 4H 


Notes. Column 1: the data cube name. Columns 2-3: the inclination angle and position angle for each ring. Columns 4—5: the mean values of the 
inclination angle and position angle shown in columns 2-3. Note that the errors only represent the standard deviations of the measured values for 
all rings. Column 6: the rotation velocity in the flat part. Column 7: the maximum rotation velocity. Column 8: the velocity dispersion in the flat 
part. Column 9: the median velocity dispersion. Columns 10-11: the ratios between the rotation velocity and the velocity dispersion. 
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Fig. 12. Best-fit of the decomposition of the circular rotation curve traced by the [Cm] line when allowing the black hole component to be free 
(left), fixed to 1.8 x 10? Mo (Feruglio et al. 2018; middle) and none (right). Black points with error bars are the circular velocities, these are the 
rotation velocities corrected for the asymmetric drift (caused by gas random motions). The solid yellow, green, brown and blue lines represent the 
black hole, stellar, gas and dark matter components, respectively. The dashed red lines are the sum of these four components. The blue diamonds 
with error bars are the circular velocities measured from the CO (9-8) line. These two lines trace the identical gravitational potential within the 


errors. 


in the quasar-host system of J2310+1855. We model the spa- 
tial distribution and extent (e.g., quantified by the half-light radii 
, and the Sérsic index) of the [C n] and the CO (9-8) lines and the 
dust continuum emission, using a 2D elliptical Sérsic function in 
Sect. 4.1. 


The extended FIR dust continuum and the CO high-J line 
emission follow a Sérsic distribution with the Sérsic index njsyq 
a bit larger than 1. Exponential dust-continuum profiles are also 
observed for nearby (e.g., Haas et al. 1998; Bianchi & Xilouris 
2011) and high-redshift galaxies (Hodge et al. 2016; Calistro 
Rivera et al. 2018; Wang et al. 2019; Novak et al. 2020). How- 
ever, the [Cu] emission profile is slightly flatter, with a Sérsic 
index of 0.59. 


Walter et al. (2022) used simple RADEX (van der Tak et al. 
2007) modeling on a z ~ 7 quasar host galaxy, and find that 
the [Cu] surface brightness in the very central part (i.e., «0.4 
kpc) is significantly suppressed by the high dust opacity, which 
increases the FIR background radiation field on the [Cn] line, 
thereby making the [C ri] fainter than one would observe with- 
out the FIR background. As shown in Fig. 9, both the dust mass 
surface density (Zaust mass > 10* Mo/kpc?) and the [C ir] underly- 
ing dust optical depth (t ~ 0.5) are very high inside the central 1 
kpc region in our quasar host galaxy. In order to test the high dust 
opacity effect on the spatial distribution of [C n], following Wal- 
ter et al. (2022), we model the [C n] line intensities as a function 
of radius with and without the FIR background determined using 
the dust emission in the host galaxy (i.e., a gray-body radiation 


of By(Taus)(1 — exp ™)). For both situations, we also consider the 
CMB as a background (i.e., a black-body radiation of B,(Tcmp); 
Tems = 2.73 x (1 +z) K). It should be noted that we only model 
collisional excitation with molecular hydrogen H2, and ignore 
collisions with electrons and neutral hydrogen Hı. We model 
the [Cn] line intensities for the whole galaxy by dividing the 
galaxy disk into nine concentric rings with ring width of 005. 
The number density of Hp is fixed to be 10?cm^?. The column 
density of [C r1] for each ring can be calculated with a fixed gas- 
to-dust mass ratio of ~10 (from our measured total dust mass 
in Sect. 4.2 and total gas mass in Sect. 4.4) and a fixed [Cu] 
abundance relative to hydrogen (7 x 107 is the best value in our 
fiting). We further assume that the average kinetic temperature 
is equal to the dust temperature at the central radius of each ring, 
and the dust temperature can be predicted from the Taust-r re- 
lation measured in Sect. 4.2. The average optical depth in each 
ring is associated with the dust mass surface density, which can 
be predicted from the Xgust_mass—/ relation measured in Sect. 4.2. 
In order to match the observed surface brightness of the [C u] 
line, we reduce the scale length of the Xau mas,-7 relation by 
a factor of 2 for the inner 0.5 kpc region, which is consistent 
with the observed surface density profiles of the dust mass for 
some nearby galaxies (e.g., Casasola et al. 2017). The full width 
at half maximum (FWHM) of the [Cn] line emission for each 
ring is measured from our observed [C n] line spectrum, which 
is extracted from each circular annulus. 
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Table 4. Derived parameters from the [C 1] dynamical modeling. 


Parameter (a) (b) (c) 
Mgg (109 Mo) (D. 29793 aa] 8 - 
Morar (10? Mc) Q) 1165 576535 — 63055 
Re: star (kpc) (D. 1.31227 idco 0219 25 
Mstar (4) 4.887347 6,661277 p RI 
Mpm (10? M5) (5) dto X500 Gd 
aco [Mo/(Kkm/spe*)] (6 0.700 0.37700 0.34790. 
Mg (1019 Mo) (T) 20595 2,021957 1.8400 
Mya (109 M5) (8) Xam 2.74228 2.65104 
p" (9) 095755 0.77208 — 0745 
Vooo (km/s) (0. Jr UID 1014116 
R»oo (kpc) (11) 1595 7. 112595" dudit 


Notes. Columns (2), (b) and (c) are for the dynamical models when allowing the black hole component to be free, fixed and none, respectively. 
» Rows 1-6: fitted parameters from the dynamical modeling in Sect. 4.4. They are the black hole mass, the stellar mass, the effective radius of the 
. stellar component, the Sérsic index of the stellar component, the dark matter mass, and the CO-to-gas mass conversion factor, respectively. Row 7: 
the gas mass (Maa, = Heol con s). Row 8: the baryonic mass (Mbary = Mstar + Meas). Row 9: gas mass fraction (feas = Meas/Mbary). Rows 10-11: 


the virial velocity and virial radius of the dark matter halo. 


» «9 the black hole mass is fixed to the value of 1.8 x 10? Mo (Feruglio et al. 2018). 


The RADEX modeling results are shown in Fig. 13. The high 
FIR background associated with the dust radiation, which has 
dust optical depth ~0.5 inside 1 kpc, indeed lowers the [C rr] sur- 
face brightness, thus changing the observed spatial distribution 
of [C u], reducing the [Cu] equivalent width, and lowering the 


^ surface brightness ratio between [C ri] and FIR emission. The ef- 


^ fect becomes weaker in the outer parts of the galaxy, as the dust 
optical depth decreases. As shown in the right panel of Fig. 9, 
the CO (9-8) underlying dust continuum optical depths are be- 
. low 0.2. And importantly, the CO (9-8) emits at longer wave- 
length than [Cn], further away from the dust emission peak. 
Thus, the FIR background radiation is not strong enough to sig- 
nificantly modify the observed CO (9-8) line intensity. In sum- 
mary, the difference of the observed surface brightness distribu- 
tion between [C rr] (with a Sérsic index of 0.59) and CO (9-8) 
(with a Sérsic index of 1.21 or 2.01) may be just the effect of 
high dust opacity (see also Riechers et al. 2014), which dimin- 
ishes the [C n] emission much more strongly in the center than 
in the outer region. 


We next compare the rotation angles (column 8 in Table 2) 
of the major axes of the tilted brightness distribution of the ISM 
tracers with the position angles (column 5 in Table 3) from the 
kinematic modeling. The photometric major axis of the [C n] line 
is close to that of its underlying dust continuum, but is almost 
perpendicular to that of the CO (9-8) line. The photometric ma- 
jor axis of the CO (9-8) line is coincident with the kinematic ma- 
jor axis of the [C n]/CO (9-8) shown in Sect. 4.3, but makes an 
angle of ~30° with respect to that of its underlying dust contin- 
uum. Such a misalignment between kinematic and photometric 
position angles is also reported for the hosts of Palomar-Green 
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quasars by Molina et al. (2021), which may indicate kinematic 
twisting. We compare the effective radii (in boldface in columns 
10-11 of Table 2) of the emission near the kinematic major axis, 
as they are less influenced by projection effects and can be inter- 
preted as the matter distribution along the kinematic major axis. 
The half-light radius of the [C n] underlying dust continuum at 
observed frame 262 GHz is larger by a factor of 1.1 than that of 
the CO (9-8) underlying dust continuum at observed frame 147 
GHz. The half-light radii of the dust continuum at both frequen- 
cies are smaller than those of their nearby lines, with ratios of 
~0.7 and ~0.5 of [C n] underlying dust continuum to the [C ri] 
line, and CO (9-8) underlying dust continuum to the CO (9-8) 
line, respectively. The half-light radius of the [C n] line is smaller 
than that of the CO (9-8) line. The more concentrated dust con- 
tinuum than the [Cn] and the CO (9-8) emission is also found 
in other z ~ 6 quasars (e.g., Wang et al. 2019; Novak et al. 
2020), high-redshift galaxies (e.g., Riechers et al. 2013, 2014; 
Chen et al. 2017; Gullberg et al. 2019), and numerical simula- 
tions (e.g., Cochrane et al. 2019; Popping et al. 2022). 


The higher compactness for the continuum than the lines re- 
flects the temperature dependence on radius; the dust surface 
brightness depends on both column density and temperature (Eq. 
1), while the lines mainly scale with temperature (e.g., Calistro 
Rivera et al. 2018). The dust continuum sizes at the rest frame 
wavelengths of 158 and 290 um are identical within ~10%, 
which is consistent with simulations conducted by Popping et al. 
(2022), who found that the dust continuum sizes remain constant 
within ~20% at observed frame wavelengths from 500 um to 
2 mm for z ~1—5 main-sequence galaxies. As these two wave- 
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Fig. 13. RADEX modeling of the [C n] line emission with (red lines) and without (blue lines) FIR background caused by high dust opacity in the 
quasar host galaxy, and compared with the observed [Cm] line emission (black circles with error bars): the [C n] surface brightness profile (left 
panel), the [C n] line equivalent width (middle panel), and the ratio between Lc, and Xp (right panel) as a function of radius. The dashed black 
line in the right panel represents the radial distribution of the optical depth at the frequency of the [Cm] underlying dust continuum. 


lengths are very close, they are both dominated by dust at similar 
temperatures. 

Our RADEX modeling of the effect of the dust opacity sug- 
gest that the actual half-light radius of the [Cu] gas is much 
smaller than that measured from the observed n — 0.59 Sérsic 
distribution listed in Table 2. The intensity of the CO (9-8) line 
depends strongly on the ISM density that drives the excitation 
of high-J CO lines. The [C n] line is mainly collisionally excited 
. by electron and hydrogen, which are heated by UV photons. The 
different effective radii of the [C n] and CO (9-8) lines may be a 
reflection of the radial dependence of the gas excitation and the 


5 radiation field. 


5.1.2. The nuclear and extended components of the ISM 


The long-wavelength dust continuum emission appears to trace 
the bulk of the star formation activity in quasar-starburst sys- 
tems at z ~ 6 (e.g., Venemans et al. 2018; Li et al. 2020c). 
Detailed SED decomposition from rest-frame UV to FIR allows 


! one to constrain the relative contribution to the continuum of dif- 


ferent components (i.e., accretion disk, torus, host galaxy; e.g., 
Leipski et al. 2013, 2014; Shao et al. 2019). We find (Fig. 8) 
that the AGN dust torus contributes ~40% of the integrated FIR 
) luminosity of J2310+1855, and the star-formation heated dust 
continuum emission dominates the FIR luminosity above rest- 
frame 50 um (Shao et al. 2019; Sect. 4.2 in this paper). The dust 
continuum decomposition results presented in Table 2 and Figs. 
B.2-B.3 reveal both a compact central component and an ex- 
tended component, which may stem from the AGN dust torus 
and the quasar host galaxy associated with star formation, re- 
spectively. The AGN dusty molecular torus is a key element in 
the AGN unification model, which obscures the accretion disk 
and the broad line region in type-2 AGN (e.g., Hónig 2019). The 
AGN torus might be unstable, warped and unaligned with the 
host galaxy orientation, and has been revealed by observations of 
several molecular lines and dust continuum in low-z AGN (i.e., 
10 pc scale; e.g., García-Burillo et al. 2014, 2021; Combes et al. 
2019). Our image decomposition assigns ~5% of the dust con- 
tinuum near both [C n] and CO (9-8) to a point source, which we 
interpret as emission from the AGN torus. 

The diameters of the dusty molecular tori measured by high- 
resolution ALMA observations of nearby Seyfert galaxies range 
from ~25 pc to ~130 pc, with a median value of ~42 pc (García- 
Burillo et al. 2021). We here consider a radius of 25 pc of the 


torus in our quasar, and with the point emission at the two wave- 
lengths listed in Table 2 and Eqs. 1 and 2 (6 fixed to be 1.90), 
we calculate a torus dust temperature of 1600 500 K and a 
torus dust mass of (1.46 + 1.14) x106M,. Esparza-Arredondo 
et al. (2021) explored the torus GDR for a sample of 36 nearby 
AGN with NuSTAR and Spitzer spectra, and found that it lies be- 
tween 1 and 100 times the Galactic ratio (i.e., 100; Bohlin et al. 
1978). From our dynamical modeling, we get an average GDR 
of ~10 for the quasar host galaxy. We assume a GDR value of 
100 for the torus alone, as this is consistent with the existence of 
gas located in the dust-free inner region of the torus. Thus, the 
molecular gas mass derived from the dust emission in the torus 
is (1.46 + 1.14) x108Mo. 

We present the best-fit CO (9—8) image decomposition re- 
sults in Fig. 7 and Table 2 with both a nuclear and extended 
components, representing the CO (9-8) emission from the AGN 
molecular torus and the AGN host galaxy, respectively. The nu- 
clear component comprises ~11% of the whole CO (9-8) emis- 
sion. The nuclear CO (9-8) emission can be converted to a 
molecular gas mass with a flux ratio between CO (9-8) and CO 
(2-1) lines (i.e., 48; Shao et al. 2019) and a CO-to-H5 conver- 
sion factor - aco, assuming Loos x Loin (Carilli & Wal- 
ter 2013). We derive an average oco value of ~0.37 from our 
dynamical study, but this may be smaller in the core. For ex- 
ample the conversion factor in the circumnuclear disk can be 
3-10 lower than that of the ISM (i.e., Usero et al. 2004; García- 
Burillo et al. 2014). Assuming a aco that is 10 times smaller 
than the overall value, the molecular gas mass is derived to be 
(2.44 + 0.28) x105 Ms, which is roughly consistent with that 
measured from the dust emission from the AGN dusty molecular 
torus. However, the decomposed unresolved emission from the 
AGN dusty molecular torus is quite uncertain given our current 
angular resolution and sensitivity. And the estimate of the molec- 
ular gas mass in the dusty molecular torus depends on many un- 
certain quantities: the torus size, the dust emissivity index, the 
GDR, the CO line ratio and the oco in the unresolved central 
region. 


5.2. The surface density of the gas and the star formation 


Next, we investigated the surface density of the gas and star for- 
mation. We measured the surface densities of gas and dust as a 
function of distance from the quasar center from the intensity 
maps using circular rings. These rings have widths of 005 and 
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0"1 for the [C n] and CO (9-8) lines, respectively, about half the 
restoring beam sizes. For each ring, we used Eqs. 1 and 2 with 
a fixed dust emissivity index of 1.90 (from the UV to FIR SED 
decomposition in Sect. 4.2) and with the dust temperature and 
dust mass surface density following the relations from our dust 
diagnostics in Sect. 4.2. Then we measured the FIR luminosity 
by integrating from 42.5 to 122.5 um and the IR luminosity by 
integrating from 8 to 1000 um. 


5.2.1. Xic] /ŻFIR 


Local spiral galaxies have a typical [C n]-FIR luminosity ratio 
of ~5 x 10? when Egg < 10!! L;/kpc?. However, the ratio is 
substantially smaller at higher FIR luminosity surface densities 
(e.g., Diaz-Santos et al. 2013, 2017). This so-called [C n]-FIR 
deficit is shown in the left panel of Fig. 14, together with re- 
sults from other spatially resolved z ~ 6 quasars (e.g., Wang 
et al. 2019) and z ~ 3 submillimeter galaxies (SMGs; e.g., Ry- 
bak et al. 2019). As shown in the right panel of Fig. 13, the ratio 
of Xjc,; and Xp of our target increases with increasing radii, 
appearing to asymptote at larger radii. Our [C n]-FIR ratios are a 
few times higher than those of spatially resolved star-forming as- 
sociations at similar FIR surface brightness. Our measurements 


») push to higher Xp than do the comparison samples, but do not 


probe deeper into the faint end (1.e., Egg < 10/0 Lo/kpc?) oc- 


. * cupied by Lyman break galaxies (LBGs; e.g., Hashimoto et al. 


' 20192) and the local luminous infrared galaxies in the GOALS 
sample (Díaz-Santos et al. 2013, 2017). 


The origin of the overall [C n]-FIR deficit of galaxies has 
. been investigated in detail (e.g., Luhman et al. 1998; Stacey et al. 


i - 2010; Diaz-Santos et al. 2013; Narayanan & Krumholz 2017; 


Lagache et al. 2018), yet remains under debate. The [Cn] line 
is mainly excited through collisions with electrons, neutral hy- 


2 - drogen Hı, and molecular hydrogen H2. Close to the quasar host 


galaxy center (corresponding to the region with the highest Xy), 
the UV photons are absorbed by large dust grains (which may be 
, reemitted at FIR wavelengths, thus boosting the FIR luminosity). 
» Therefore, the hydrogen and electrons will be less heated, giv- 
ing rise to few collisions to excite ionized carbon. In addition, 
the dust absorption reduces the number of ionizing photons, de- 


' creasing the extent to which carbon becomes ionized (e.g., Luh- 


man et al. 2003). Other popular explanations for the [C n]-FIR 
deficit are the saturation of [C n] emission at very high tempera- 
tures (e.g., Tas > Tic; e.g., Muñoz & Oh 2016) or in extremely 
dense environments (in which more carbon is in the form of CO 
than Ct; e.g., Narayanan & Krumholz 2017). 


The observed larger [C n]-FIR deficit in the center of our 
quasar host galaxy may be a reflection of the high dust temper- 
ature at small radius, and the different dependence of Xpy and 
Zc ON Taust: AS shown in the right panel of Fig. 14, Xy is an 
integral over the gray-body relating to the dust temperature with 


a power-law index (~ T3004. black line), larger than that (~ T^) 


dust ? 
of the Stefan-Boltzmann law. However, Xic scales as ~ T2 


(green line). Thus, higher dust temperatures in the galaxy cen- 
ter produce a larger [C n]-FIR deficit (see also Gullberg et al. 
2015; Walter et al. 2022). In addition, as demonstrated by Wal- 
ter et al. (2022) and our RADEX test in Sect. 5.1.1, the high dust 
opacity, which increases the FIR background radiation field on 
the [C un] line (see also Riechers et al. 2014), suppresses the ob- 
served [C 1] line, especially in the central region, thus enhancing 
the [C n]-FIR deficit. 
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5.2.2. XM. and XsER 


gas 


The Kennicutt-Schmidt (KS) power-law relation between the gas 
mass and the SFR surface density, Xspr ~ DR , describes how 


efficiently galaxies turn their gas into stars, which enables us to 
understand galaxy formation and evolution across cosmic time. 
The KS relation is nearly linear (N —1-1.5) in galaxies ranging 
from the local Universe to high redshifts (z « 4; e.g., Kennicutt 
1998; Bothwell et al. 2010; Genzel et al. 2010, 2013), indicat- 
ing that the star formation processes may be independent of cos- 
mic time. It is still unknown whether it holds for higher-redshift 
galaxies up to the reionization epoch. 

In order to determine whether the spatially resolved ISM of 
J2310+1855 at z = 6.00 deviates from the star formation law 
observed in local galaxies and high-z star-forming galaxies, Fig. 
15 plots the SFR surface density as a function of the molecu- 
lar gas mass surface density of the resolved ISM of J2310+1855 
(black open circles and diamonds with error bars), and of other 
populations from local to high-z galaxies in the literature (refer- 
ences in the caption). The molecular clouds where stars form 
are primarily composed of molecular hydrogen (H5) and he- 
lium. Given cosmological abundances, we consider the mass of 
the star-forming gas to be the mass of the molecular hydrogen 
times 1.36. We measure the gas mass surface density Xy... in 
each ring from the CO (9-8) intensity map with a flux ratio be- 
tween CO (9-8) and CO (2-1) lines of ~8 (Shao et al. 2019) and 


assume Loon RI E cos (Carilli & Walter 2013) and a con- 


version factor aco ~ 0.37 Mo/(K km/s pe?) from the dynamical 
modeling in Sect. 4.4. In addition, we also measure Xy, from 
the [Cr] intensity map with ajc ~ 3.2 Mo/Lo from the dy- 
namical modeling in Sect. 4.4. We similarly measure the SFR 
surface density Xsrg from the resolved dust continuum by con- 
verting the IR luminosity to SFR using the formula in Kennicutt 
(1998). Note that we subtracted the IR luminosity contributed by 
the AGN torus as a point source (see Sects. 4.1 and 5.1.2). 

Our SFR and gas surface densities are at the high end of the 
samples shown in Fig. 15, comparable to the SMGs at 2 < z < 4. 
The best-fit KS relation with Xy,,, measured from the CO (9-8) 
line for J2310+1855 is i 


log;o XsrR 


login Èm, 
Meine = —20.50(41.87) + 2.50(+0.21) x £10 2M; 
o 


—— m (9 
Mo /kpc? 9» 
and the best-fit KS relation with È Meas measured from the [C n] 
line is 


log;o XsrR 


Og10 Xu. 
Me/yr/kpc? 


l 
= —10.59(+0.64) + 1.37(+0.07 1 
0.59(+0.64) + 1.37(+0.07) x Moe (10) 


for which we do not fit with the inner most measurement, as the 
observed [C ir] emission in the central of the galaxy is highly 
suppressed due to the high dust opacity (see Sect. 5.1.1). 

The power-law index N = 2.50 + 0.21 measured from the 
CO (9-8) line is higher than that of local disk galaxies (e.g., 
Kennicutt 1998; Kennicutt et al. 2007) and high-redshift star- 
forming populations (e.g., Genzel et al. 2013; Sharon et al. 
2019). This may indicate that this quasar host galaxy at z = 6.00 
is undergoing much faster star formation than do galaxies in 
the local Universe and at lower redshifts. However, we should 
note that, as mentioned in other studies (e.g., Daddi et al. 2010; 
Thomson et al. 2015) the CO-to-H» conversion factor aco can 
vary for different galaxy types, position scales, and metallici- 
ties. We adopted the oco value from our dynamical studies in 
Sect. 4.4, which only represents the global value, and ignored 
any possible variation within the galaxy. This will bias the slope 
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Fig. 14. [C u]-FIR deficit. Left panel: [C n]/FIR luminosity ratio as a function of FIR surface brightness. The red open diamonds with error bars 
connected by the solid red line represent measurements for J2310+1855. The comparison samples are spatially resolved z ~ 6 quasars (green 
triangles connected by dashed lines; Wang et al. 2019) and z ~ 3 SMGs (purple triangles connected by dashed lines; Rybak et al. 2019), other 
z > 6 quasars (yellow stars; Wang et al. 2013; Venemans et al. 2016, 2017; Decarli et al. 2018), other high-z SMGs (pink triangles; Riechers et al. 
2013; Neri et al. 2014; Gullberg et al. 2018; Rybak et al. 2019), the GOALS sample of local luminous infrared galaxies (gray circles; Díaz-Santos 
et al. 2013, 2017), and LBGs (green triangles; Capak et al. 2015; Jones et al. 2017; Hashimoto et al. 20192). Right panel: Comparison between 


. ofthe KS law we found in Fig. 15. In addition, the indirect mea- 
surements of the low-J CO line transitions requires assuming 
a CO excitation model, which brings additional uncertainties 


( > when comparing KS relations from different studies. We sim- 
*$ ply use a bulk ratio between the ALMA CO (9-8) emission and 
I the VLA CO (2-1) emission, which in reality will decrease with 


* radius. This leads to an over-estimation of the slope of the KS 
law. Low-J CO observations with higher spatial resolution will 


== be needed to improve on this analysis. And there appears to be 


an excess in CO (9-8) luminosities in distant starbursts due to 
for example shock excitation (e.g., Riechers et al. 2021a), which 


4 - would make the CO (9-8) line a poorer tracer of gas mass. 


The power-law index N = 1.37 + 0.07 measured from the 
= [Cu] line is consistent with both the local and other high-redshift 
samples (e.g., Bothwell et al. 2010; Carilli et al. 2010; Genzel 


FI et al. 2010; Thomson et al. 2015). This may indicate that this 


quasar host galaxy at z — 6 in our work follows similar star 
formation process with that of the local Universe and lower red- 
shifts. However, the value of the o(c,; might be variable within 
the galaxy. The depletion times defined by 2y,,,/2srr for the 
ISM of our quasar host galaxy are 1-100 Myr, increasing from 
the inner region (corresponding to the highest Xy,,, and Xsrg) 
to outer. The short gas consumption timescales suggest a rapid 
starburst mode for our quasar host galaxy. 


5.3. The asymmetric line profile 


The integrated spectra of the gaseous disks of many galaxies 
both in the local Universe and at high redshift often have a char- 
acteristic double-horned shape (e.g., Roberts 1978; Walter et al. 
2008; Shao et al. 2017). Higher velocity dispersion and steeper 
density profiles of the gas tend to decrease the depth of the val- 
ley between the horns, and can even reduce the global profile to 
flat-topped Gaussian shapes (e.g., de Blok & Walter 2014; Stew- 
art et al. 2014). In addition, the gaseous disk can be disturbed by 
interactions of galaxies with their neighbors and environments 


jj the trends of the surface brightnesses of [C n] (green circles with error bars) and the scaled FIR (black circles with error bars) luminosities, as a 
» function of dust temperature. The green and black lines are the power-law fits to the Xjcm and Xp as a function of Taust, respectively. 


(e.g., van Eymeren et al. 2011; Bok et al. 2019; Watts et al. 2020; 
Reynolds et al. 2020), ram-pressure stripping (e.g., Gunn & Gott 
1972; Kenney et al. 2015) and gas accretion from the cosmic 
web (e.g., Bournaud et al. 2005). These distortions may drive 
morphological (e.g., the gas distribution) and kinematic (e.g., the 
velocity field, and the global spectral profile) asymmetries. We 
now turn to the spectral profile of J2310+1855 to understand it 
in the context of these ideas. 


The [Cn] and CO (9-8) spectral profiles of J2310+1855 are 
asymmetric as shown in Figs. | and 2. The spectra are slightly 
enhanced on the red side. The enhancements ([S) redpeak — 
S y bluepeak]/S'y,bluepeak3 from the double-Gaussian fitting) are 
about 26% + 10% and 18% +14% of the [C n] and CO (9-8) lines, 
respectively. The two possible peaks are separated by 213 + 9 
and 215 + 14 km/s for the [C n] and CO (9-8) lines, respectively. 
From our kinematic modeling of the [Cn] and CO (9-8) data 
cube in Sect. 4.3, we found a median gas dispersion of 54*6 and 
645 km/s for the [C n] and CO (9-8) lines, respectively (Table 
3). The ratios between the separation of the two peaks in the red 
and blue parts and the median velocity dispersion are 39 and 
39707 for the [C n] and CO (9-8) lines, respectively. In addition, 
as shown in Sect. 4.3, the radial distribution of [C ir] and CO (9- 
8) lines are not steep (e.g., the Sérsic index is around 1). Thus, 
we are able to detect the double-horned profiles for both the [C ri] 
and CO (9-8) lines (we consider the two peaks to be resolved if 
the separation of the two peaks is twice the velocity dispersion). 
The redward enhancement of the spectra may be due to more 
material on one side of the galaxy, or temperature differentials 
from one side to another. 


5.4. lonized and molecular gas kinematics 


We studied the gas kinematics with both [C n] and CO (9-8) lines 
using "PBARoro in Sect. 4.3. These two lines present consistent 
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| Fig. 15. Spatially resolved KS plot for J2310+1855 (black open circles with error bars, and its best-fit power-law shown as black solid line 
when calculating Xy, from the CO (9-8) line; black open diamonds with error bars, and its best-fit power-law shown as dashed black line when 
calculating Zm, from the [Cn] line), compared with the surface density KS relations of other galaxy samples including local starburst galaxies 
(brown filled circles; Kennicutt 1998), local spiral galaxies (blue filled circles; Kennicutt 1998; Bigiel et al. 2010), z ~1-3 star-forming galaxies 
(pink filled circles; Genzel et al. 2010; Freundlich et al. 2013; Tacconi et al. 2013), z ~2-4 SMGs (blue filled squares; Bothwell et al. 2010; Genzel 
~ et al. 2010; Carilli et al. 2010), star-forming dusty galaxies (blue open stars; Villanueva et al. 2017), local blue compact dwarf galaxies (red open 
triangles; Amorín et al. 2016), and some spatially resolved galaxies — Milky Way clumps (green open diamonds; Heiderman et al. 2010; Evans et al. 
2014), star-forming regions of the spiral galaxy NGC 5194 (gray filled circles; Kennicutt et al. 2007), the lensed SMG SDP.81 at z = 3.042 (red 
pluses; Hatsukade et al. 2015), the young low-metallicity starburst galaxy RCSGA 032727 at z = 1.7 (brown crosses; González-López et al. 2017), 
massive star-forming galaxies around z = 1.2 (purple open diamonds; Freundlich et al. 2013), the z = 1.5 star-forming galaxy EGS13011166 
(dashed red line; Genzel et al. 2013), the lensed SMG SMMJ21352 at z = 2.3 (green filled triangles; Thomson et al. 2015), and the H u regions 
in the nearby spiral galaxy M33 (black crosses; Miura et al. 2014). The solid gray lines represent gas depletion timescales of 1Myr, 10 Myr, 100 
Myr, and 1 Gyr. We applied a 1.36 correction factor to the molecular hydrogen mass for the presence of helium for the comparison samples. For 
our target, the conversion factor to the gas mass is from the dynamical modeling, which should include the contributions from all other elements 
(e.g., helium). When deriving the index of the KS power-law with the X Meas derived from the [C ir] line, we do not use the innermost measurement, 

as the observed [C n] emission in the central of the galaxy is highly suppressed due to the high dust opacity (see Sect. 5.1.1). 


kinematic properties, including the gas disk geometry and the for the two lines. The modest inclination angle (~40°) makes 

gas motions. J2310+1855 an ideal target to investigate the gas kinematics 
in the early Universe. The rotation speed as well as the veloc- 
ity dispersion (discussed below) are also consistent between the 

5.4.1. [Cu] versus CO (9-8) [C ir] and CO (9-8) lines. This suggests that the [C i] and CO (9- 
8) gas are coplanar and corotating in the host galaxy of quasar 

The [Cu] and CO (9-8) emission trace similar gas disk geome- J2310+1855. 

tries in our quasar J2310+1855 host galaxy. As shown in Table The rotation speed of the [C n] line presented in the left panel 

3, the inferred inclination angle and position angle are identical of Fig. 11 is roughly constant (~250 km/s) with radius. How- 
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ever, the rotation speed of the CO (9-8) line at ~1.2 kpc is some- 
what lower than that of the other two rings. The measured rota- 
tion speed of the [C ri] line is overall consistent with that of the 
CO (9—8) line, given that the errors in Fig. 11 do not include the 
covariance with other fitting parameters. The [C n] and CO (9-8) 
data cubes used in the kinematic modeling in Sect. 4.3 are from 
the same telescope and follow an identical data reduction pro- 
cess. They have a similar velocity resolution (^65 km/s), which 
is sufficient to sample the intrinsic width of both lines, and the 
line sensitivities are also comparable (i.e., S/Npeak > 10). The 
only difference is the angular resolution: that of the [Cn] data 
is roughly two times better than that of the CO (9-8) data. Con- 
sidering that we measured the gas rotation speed in rings with 
widths that are half the angular resolution for both lines, and 
that "PBAROoLo performs well even when the galaxy is resolved 
with a small number of beams (Di Teodoro & Fraternali 2015), 
the spatial resolution difference may not affect the inferred ro- 
tation curve significantly. Thus, the lower rotational speed of 
the ~1.2 kpc ring from CO (9-8) kinematical modeling may be 
due to higher random motions of the CO (9-8) gas in that ring. 
The maximum rotation velocities (~250 km/s) of J2310+1855 
in this work are smaller than that (~400 km/s) of quasar ULAS 
J131911.29+095051.4 at z = 6.13 in our previous work (Shao 
et al. 2017), but are at the high end of the rotation velocities of a 
sample of 27 z > 6 quasars (Neeleman et al. 2021). The velocity 

ş dispersions of the [Cm] and CO (9-8) lines have consistent me- 
ag dian values of 54*6 and 6475 km/s, respectively, as shown in the 

» middle panel of Fig. 11 and Table 3. This is different from what 
, is normally observed in nearby galaxies where the ionized gas 
(i.e., Ha) velocity dispersions are substantially higher than the 
molecular gas (i.e., CO) velocity dispersions (e.g., Fukui et al. 


J 2009; Epinat et al. 2010). However, because the angular resolu- 
YJ tion of these two lines are different, the velocity dispersions are 


measured over different volumes, and in the Milky Way, the ve- 
locity dispersion in molecular clouds is proportional to the size 
and mass of the clouds (e.g., Heyer & Dame 2015). In addition, 
our velocity dispersions are much smaller than the average ve- 
locity dispersion (i.e., 129 + 10 km/s) estimated using ALMA 
[C u] data at a resolution of ~0’’25 in a sample of 27 z > 6 undis- 


f > turbed quasar host galaxies (Neeleman et al. 2021). 


There is a general trend that the gas velocity dispersion in- 
creases with redshift or cosmic time (e.g., Glazebrook 2013). 
The predicted velocity dispersion for z = 6 galaxies is about 
| ^80 km/s following the empirical relationship for 0 < z < 4 
galaxies with measured velocity dispersions from either ionized 
or molecular gas (Übler et al. 2019). The predicted velocity dis- 


persion is consistent with our measured ones — 36*};-8078 km/s 


for the [C n] line and 54*3-98*6 km/s for the CO (9-8) line. This 
may indicate that the luminous AGN in the center has little in- 
fluence on the velocity dispersion of the quasar host galaxy. The 
physical process to drive and maintain the velocity dispersion 
over cosmic time is still not clear. But a constant energy input 
is required to retain the turbulence in the ISM, or it will decay 
within a few megayears as proposed by theoretical works (e.g., 
Stone et al. 1998; Mac Low et al. 1998). This energy supply may 
be associated with either stellar feedback (i.e., winds; expand- 
ing H u regions; supernovae) or other instability processes (i.e., 
clump formation; radial flows; accretion; galaxy interactions). 
The gas outflow along the line of sight traced by the blueshifted 
absorption of the OH* line would be very much in line with help- 
ing the upkeep of turbulence. 

The rotation-to-dispersion ratio Vio;/c: is a measure of the 
kinematic support from ordered rotation (i.e., circular motions) 
versus random motions (i.e., noncircular; turbulence) in a sys- 


tem. We make two kinds of measurements of V,.¢/o in Table 3 
— the ratio between the rotation velocity and velocity dispersion 
in the flat part of the rotation curve (column 10) and the ratio 
between the maximum rotation velocity and the median velocity 
dispersion (column 11). The ratios are above 4 for both the [C u] 
and CO (9-8) lines, which are higher than the cutoffs of a rotat- 
ing system (i.e., a cutoff of 1 from Forster Schreiber et al. 2009; 
Epinat et al. 2009; a cutoff of 3 from Burkert et al. 2010; Fórster 
Schreiber et al. 2018), indicating that the [C ir] and CO (9-8) gas 
motion is dominated by rotation. These V,o/o values are typical 
of other z ~ 6 quasar host galaxies traced by resolved [C n] in 
Neeleman et al. (2021). Similarly, our measured V,ot/0 ratios are 
comparable to those of z — 2 star-forming galaxies (e.g., Gen- 
zel et al. 2008; Fórster Schreiber et al. 2009, 2018; Wisnioski 
et al. 2015), and z ~ 4.5 gravitational lensed dusty star-forming 
galaxies (e.g., Rizzo et al. 2021). 

The inclination angle measured from our kinematic model- 
ing is larger than the one (25?) derived by Tripodi et al. (2022), 
who also used ?PBARoro. They constrain the inclination angle 
«30? considering that Mayn should be larger than My,, which is 
derived from the CO (2-1) emission line (Shao et al. 2019) and 
adopting an o, value of 0.8 Mo /(K km/s pc’). As noted in Sect. 
4.4, the values of a. may differ from source to source (e.g., 
Mashian et al. 2015), and based on our dynamical model, we de- 
rive an a, value of 0.37 for J2310+1855. The loose constraint 
ON Qo (i.e., we adopt a wide value range of 0.2-14; Mashian 
et al. 2015) allows for a more flexible fitting of the inclination 
angle in the PBAnoro modeling. The initial guess of the incli- 
nation angle of 40° in the ?PBAnoro modeling is from the ratio 
between the photometric minor and major source sizes (which 
are roughly along the kinematic minor and major axes, respec- 
tively) of the CO (9-8) emission assuming an intrinsic round gas 
disk (before inclination). However, the observed kinematic mi- 
nor source size is larger than the kinematic major source size 
(after inclination) for the [C n] line, which indicates that the in- 
trinsic kinematic minor source size of [C n] must be larger than 
its intrinsic kinematic major source size (before inclination). As 
we discuss above, the [C n] and CO (9-8) gas are coplanar and 
corotating in the host galaxy of quasar J2310+1855. Thus, [C 11] 
and CO (9-8) gas are well mixed. Therefore, the intrinsic kine- 
matic minor source size of CO (9-8) might be larger than its in- 
trinsic kinematic major source size. As a result, from the kine- 
matical aspect, the inclination angle should be larger than 40^, or 
conservatively the inclination angle cannot be much below 40° 
considering a 5% uncertainty of the initial guess of the inclina- 
tion angle. 


5.4.2. Outflow traced by OH* 


The OH* line has been detected in emission, absorption and 
P-Cygni-shaped profiles in galaxies (e.g., van der Werf et al. 
2010; Rangwala et al. 2011; Spinoglio et al. 2012; González- 
Alfonso et al. 2013; Riechers et al. 2013, 2021a,b; Berta et al. 
2021; Butler et al. 2021; Stanley et al. 2021). The production 
of OH* is mainly driven by H* + O — O* +H followed by 
O* + H? — OH* +H and H* + OH  OH* +H (e.g., van der 
Werf et al. 2010). The OH* is in a unstable state, and it reacts 
rapidly with Hz molecules following the main chemical reac- 
tions: OH+ + Hy — H50* +H and then H20* +H; —^ H30* +H 
(e.g., Gerin et al. 2010). The key initiation for this simple chem- 
ical sequence is the production of H*, and it is argued that the 
atomic H ionization is most likely dominated by X-rays and/or 
cosmic rays (e.g., González-Alfonso et al. 2013). The OH* emis- 
sion shows double-horned profile with comparable intensities of 
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^ 0.30 + 0.05 mJy of the two peaks. Our measured luminosities 
of the OH* emission and FIR follow the weak correlation be- 
tween these two parameters found by Riechers et al. (2021a). At 
—384 + 128 km/s (relative to the [C rr] redshift from our ALMA 
Cycle 0 observations; Wang et al. 2013), we detect a ~3o dip 
with a peak value of —0.18 + 0.05 mJy. This absorption velocity 
shift is at the high end of the OH* absorption velocity shifts in 
the range of ~130-360 km/s of z=2-6 starburst galaxies (Riech- 
ers et al. 2021a), and is comparable to the widths of the [C ri] 
and CO (9-8) lines of J2310+1855. This P-Cygni profile (i.e., 
blueshifted absorption and redshifted emission) may indicate the 
outflow of gas along the line of sight. The outflow signature of 
broad blue-shifted absorption line of J2310+1855 is recently re- 
ported by Bischetti et al. (2022) through VLT/X-shooter obser- 
vations. The OH* absorption and emission components are co- 
spatial (left panel of Fig. 3). This is different from, for example, 
a hyperluminous dusty starbursting major merger ADFS-27 in 
Riechers et al. (2021b), where the offset is >1 kpc. 

In order to further understand the outflow traced by OH* ab- 
sorption, we estimate the outflow mass with the observed OH* 
absorption spectrum and its best-fit Gaussian model. The opti- 
cal depth of an unsaturated absorption line can be calculated 
with Tine = —1N(firans), Where fuas, is the fraction of the con- 
tinuum emission that is still transmitted, which can be expressed 
as (S von +S vine) /S veon Where S yin is the flux density of the con- 
tinuum subtracted OH* absorption, and Sy, is the flux density 
of the continuum emission for which we use a constant value 
of 1.53 mJy for the whole line. We measure the peak optical 
depth rog- = 0.13 + 0.04. The velocity-integrated optical depths 
Tou dv are 35 + 7 and 44 + 31 km/s by integrating the OH* ab- 
sorption region and the best-fit Gaussian region of the OH* ab- 
sorption, respectively. The OH* column density Non: is related 
to the velocity-integrated optical depth rog. dv through rog dv = 
APMANog- /8/n, where A is the Einstein coefficient with a value of 
2.11 x 107? s7! (e.g., Butler et al. 2021), and A is the rest frame 
wavelength of the OH* line (i.e., 0.029 cm). We thus derive Nop- 
to be (1.7 x 0.3) x 10!^ and (2.2 + 1.5) x 10'^ cm”, respectively. 
Then assuming a relative abundance of log,,(Nou+/Nu) = —7.8 
(e.g., Bialy et al. 2019), we measure the neutral hydrogen col- 
: umn density Ng of (1.1 + 0.2) x 10” and (1.4 + 0.9) x 10? 
cm^?, respectively. With a source radius of 1.3 kpc (the effective 
radius of CO (9-8) line in Table 2), we thus derive the neutral 
hydrogen mass and neutral gas mass (i.e., with a correction fac- 
tor of 1.36 for the presence of helium) to be (4.6 + 0.9) x 108 
and (6.2 + 1.2) x 108, and (5.8 + 4.0) x 10? and (7.9 + 5.5) x 108 
Mo, respectively. The neutral gas mass of the outflow is 3 + 1% 
and 4 + 3% of the molecular gas mass in the host galaxy, respec- 
tively. Our neutral hydrogen mass is a few times smaller than the 
median value of ~ 2 x 10? M, for a sample of 18 z=2—6 starburst 
galaxies that are also traced by OH* absorption (Riechers et al. 
20212). 

Finally, we do not find evidence of gas outflow in the spectra 
of [C i] or CO (9-8) emission lines in the form of high-velocity 
wings (see Figs. 1 and 2). However, Tripodi et al. (2022) reported 
the presence of high-velocity components on both the red and 
blue sides of the [C n] emission line in J2310+1855. This dis- 
crepancy might be due to differences in the sensitivities and/or 
angular resolution of the data used in the analysis. We have better 
angular resolution (0711»x0"09 versus 017x015), although the 
sensitivities are comparable (0.17 mJy/beam per 17 km/s versus 
0.23 mJy/beam per 8.5 km/s). To understand this discrepancy, we 
used a nature weighting (robust=2.0) to TCLEAN our [C u] data, 
which give more weights on short baselines. The resulted clean 
beam size is 016 x 0713, which is similar to the one of Tripodi 
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et al. (2022), and a sensitivity of 0.15 mJy/beam per 17 km/s. 
The residual spectrum between the observed one and the best-fit 
double-Gaussian model is only noise-like, with an rms of 0.70 
mJy. The reality of the high-velocity wings reported by Tripodi 
et al. (2022) is therefore doubtful. As shown in previous studies, 
the [C ir] emission line has proven to be a difficult tracer of out- 
flow activity (e.g., Meyer et al. 2022) and, as shown by Spilker 
et al. (2020), strongly lensed dusty star-forming galaxies at z > 4 
that show clear OH outflow do not display high-velocity wings 
in their [C rr] emission lines. In addition, our CO (9—8) data does 
not have enough S/N to look for broad line wings. 


5.5. The dynamics of the quasar-host system 


The dynamical modeling of the rotation curves can yield de- 
tailed knowledge of the composition and distribution of mat- 
ter in the quasar-host system. In Sect. 4.4, we decomposed the 
J2310+1855 quasar-host system into four dynamical compo- 
nents (black hole, stars, gas, and dark matter) and modeled the 
potential contribution of each component to the measured circu- 
lar velocities of the high-resolution [C n] data cube. 


5.5.1. The dynamical structures 


The black hole mass is a physical parameter that is critical for 
understanding the coevolution of the central SMBHs and their 
host galaxies from the local to the high-redshift Universe (e.g., 
Kormendy & Ho 2013; Shen 2013). The mass of the black hole 
of the quasar J2310+1855 at z = 6 is on the order of 10? Mo 
(Jiang et al. 2016; Feruglio et al. 2018), which is measured us- 
ing the width of NIR emission lines and relationships calibrated 
with low-redshift AGN. The black hole component is important 
to shape the gravitational potential in the central region of the 
host, contributing ~120 km/s (for Mgy = 10? Mo) at a radius 
of 0.29 kpc (the innermost position we can reach with our cur- 
rent resolution), where we measured a circular velocity of ~250 
km/s. For our target, the value of Re. star is strongly covariant with 
the black hole mass and the inferred physical parameters of the 
stellar component. As shown in Fig. 12 and Table 4, a factor of 
2 change in the black hole mass causes a factor of 4 change in 
the inferred stellar mass. If we consider an extreme case without 
a central black hole, the stellar profile needs to become strongly 
centrally concentrated, with a stellar effective radius of 0.214933 
kpc to make up for the black hole (as shown in the right panel 
of Fig. 12 and column (c) of Table 4). There is no direct obser- 
vations to constrain the distribution of the stellar component in 
this quasar host galaxy. A highly concentrated stellar component 
might be possible, and thus we cannot definitely remove such an 
extreme case without a central SMBH. When making the black 
hole mass a free parameter, we confirm that the black hole mass 
is on the order of 10? Mo; this is the first time that a black hole 
mass has been dynamically measured at z ~ 6. The evidence of 
the existence of a central SMBH is marginal, and we need more 
and better data to verify (i.e., James Webb Space Telescope ob- 
servations to constrain the stellar property, and higher-resolution 
ALMA [Cu] observations). However, we can built a strong case 
for the existence of a > 10? Mọ black hole from our dynami- 
cal modeling, unless the stellar distribution in the case without 
a black hole is strange (but not impossible), as do the other line 
width metrics (e.g., Jiang et al. 2016; Feruglio et al. 2018). In ei- 
ther scenario (making the black hole mass free or fixed) the gas 
mass stays stable, as the effective radius of the gas component 
is fixed based on our ALMA CO (9—8) observations and the gas 
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mass is mainly determined by the circular velocities in the outer 
regions. 


The mass and size of the stellar component can provide es- 
sential insights into galaxy evolution (e.g., Hodge et al. 2016). 
The stellar mass and size correlation of galaxies across cos- 
mic time provides an important view of the assembly history 
of galaxies (e.g., Kawinwanichakij et al. 2021). However, for 
quasar host galaxies at z — 6, it is very difficult to determine 
these stellar properties through rest frame optical and NIR imag- 
ing, as the starlight is overwhelmed by the central quasar. Thus, 
the comparison with the z ~ 6 objects on the stellar mass-size 
plane is currently not feasible. This also leads to limitations 
on our dynamical modeling including a stellar component of 
the quasar-host system in Sect. 4.4. Our modeled stellar mass 
(around 10? M,) has a large uncertainty. It has been found that 
not all stellar components are necessarily co-spatial with the 
gas/dust for SMGs. Thus, the total Mstar may be under-predicted 
with the dynamical modeling based on the kinematics of the 
gas (e.g., Gómez-Guijarro et al. 2018). High-resolution stellar 
continuum emission observations toward our target and simi- 
lar luminous quasars in the early Universe by for example the 
James Webb Space Telescope (JWST) are critical to investigate 
the Mstar — Re; sir relation and the role of these z ~ 6 quasars in 
galaxy evolution. The stellar effective radius Re: star (in Table 4) 
is roughly similar to our measured CO and [C n] effective radii, 
but larger than the dust continuum effective radius, albeit with 
) large uncertainties (in Table 2). This is consistent with the simu- 
lated TNG50 star-forming galaxies (Popping et al. 2022). Figure 
12 shows dynamical evidence for a massive (^10? Mo) stellar 
component, especially when we fix the black hole mass to be 
1.8 x 10? Ms. In this case, the contribution to V, from stars rises 
. monotonically to the center of the galaxy, as is typical for nearby 
massive spiral galaxies (e.g., Simard et al. 2011; Gao et al. 2019, 
2020), and such a stellar bulge signature has already been proven 
to exist using a dynamical method in some dusty star-forming 
galaxies at z ~ 4 — 5 (e.g., Rizzo et al. 2020, 2021; Tsukui & 
Iguchi 2021). 


The gas component is well constrained by our ALMA 
CO (9-8) image decomposition and the dynamical modeling 
(see Sects. 4.1 and 4.4). The gas fraction (fzas = Meas/Mbary 
and Mya, = Mia + Meas) is >70%, and thus it dominates 
the quasar-host system. The large reservoir of molecular gas 
(Meas ~ 10!°Mo) is responsible for fueling the bulk of the star 
formation in the central region of the quasar host galaxy. We 
should point out that the scale radius of the gas component in the 
dynamical modeling is from the CO (9-8) intensity map assum- 
ing an exponential gas distribution. However, the overall molec- 
ular gas (e.g., traced by CO (1—0)) might be more extended than 
that traced by high-J CO. We adopt a 1.5 times larger gas scale 
radius to do the dynamical modeling, and get an oco value of 


0.16 
0.737035. 


The dark matter component is less constrained. The mass of 
the dark matter halo is 16+76% and 9128% of the total dynamical 
mass of the quasar-starburst system within a radius of two times 
the gas half-light radius, when making the black hole mass free 
and fixed, respectively. The virial velocity and radius of the dark 
matter halo both have large uncertainties as shown in Table 4. 
The virial radius of the dark matter halo is much more extended 
than the gas distribution as shown in Table 2, which is consistent 
with numerical simulations (e.g., Khandai et al. 2012). 


5.5.2. The coevolution of the quasar and its host galaxy 


To study the coevolution of central black holes and their host 
galaxies, we typically investigate the evolution of the Mgy — c 
and the Mgy — Mbuige relations across cosmic time (e.g., Kor- 
mendy & Ho 2013), where c and Mpulge are the velocity disper- 
sion and total mass of the bulge, respectively. 

Considering that the stars and molecular gas trace the same 
underlying potential, and following Davis et al. (2013), we as- 
sume that the second moments of the stars and molecular gas are 
equal. Thus, we can derive a simple zeroth-order estimate of the 
velocity dispersion: 


2 2 — y2 2 
V stars + Ostars = V zas + T gas? 


(1) 


where V,4, and c4, are the velocity and velocity dispersion of 
the molecular gas, for which we use the maximum circular ve- 
locity (~245 km/s; the maximum rotation velocity corrected by 
the asymmetric drift) and the median velocity dispersion (~64 
km/s) of the CO (9-8) line shown in Table 3 as proxies. Vetars is 
the stellar velocity, for which we use the modeled stellar veloc- 
ity (40 and 99 km/s when making the black hole mass free and 
fixed, respectively) at the half-mass radius of the stellar com- 
ponent (1.76 and 1.39 kpc when making the black hole mass 
free and fixed, respectively). The predicted stellar velocity dis- 
persions (stars) are 250 and 233 km/s where we free or fix the 
black hole mass, respectively. Based on the local relationship for 
nearby galaxies from Kormendy & Ho (2013), we predict a bulge 
velocity dispersion of ~335 +20 and ~300+ 16 km/s for a central 
black hole of mass 2.97 x 10? and 1.8 x 10? M, (fixed value in 
the dynamical modeling), respectively. These bulge velocity dis- 
persions from the local relationship are higher than the predicted 
stellar velocity dispersions including the scatter in the local re- 
lationship, also when we adopt the [Cu] kinematic parameters. 
This may indicate that the black hole evolves faster than the host 
galaxy of J2310+1855. 

We present the baryonic mass Mpary from the best-fit dynami- 
cal modeling and the ratio of Mgy/Mbpary between the black hole 
mass and the baryonic mass as a function of radius in Fig. 16. 
Inside a radius that is two times the half-light radius of the CO 
emission (72.6 kpc), the measured ratio of Mgy/Mbary is -28'? 
and ~16*5 times higher than the local relationship when making 
the black hole mass free and fixed in the dynamical modeling, 
respectively. If we consider a stellar-spheroidal component (i.e., 
the stellar component in our dynamical modeling in Sect. 4.4), 
the Mgu/Mgtar will be >50 times higher than the local relation- 
ship (Kormendy & Ho 2013). This may indicate that the central 
SMBH grows the bulk of its mass before the formation of most 
of the stellar mass in this quasar host galaxy in the early Uni- 
verse. Wang et al. (2013) and Tripodi et al. (2022) also found the 
same result for quasar J2310+1855. This is consistent with what 
has been found for other high-luminosity z ~ 6 quasars (e.g., 
Wang et al. 2019). 

However, as we discussed in Sect. 5.5.1, the black hole mass 
and the parameters (e.g., the mass, the effective radius, and the 
Sérsic index) of the stellar component are somewhat degenerate 
in the dynamical modeling. Thus, the modeled black hole mass, 
stellar mass and stellar velocity adopted in the discussion above 
have substantial and coupled uncertainties to the investigation of 
the quasar-host coevolution. Higher angular resolution and sen- 
sitivity observations toward the gas and stellar components are 
needed to better constrain these parameters with dynamical mod- 
eling. In addition, Eq. 11, which we used to predict the stellar 
velocity dispersion comes with additional uncertainties. Davis 
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Fig. 16. Baryonic mass Mya (red lines) measured from the dynamical 
modeling in Sect. 4.4 as a function of radius, and the corresponding 
ratio (green lines) between the black hole mass Mgg and the baryonic 
mass My. The solid and dashed lines represent the scenarios in which 
the black hole mass is free and fixed during the dynamical modeling, 
respectively. 


et al. (2013) have compared the predicted stellar velocity disper- 
sion using this equation with the measured values. Some data 
points with stellar velocity dispersion similar to the circular ve- 
locity (which is fairly common in elliptical galaxies) appear as 
obvious outliers, which may indicate that the zeroth order ap- 
proximations may break down in this situation. 


6. Summary 


We have reported on ALMA high-resolution (sub-kiloparsec- 
to kiloparsec-scale) observations of the [Cu], CO (9-8), and 
OH"* (1,01) lines and their underlying dust continuum toward 
, the FIR luminous quasar SDSS J2310+1855 at z = 6.0031. We 
} have investigated the ISM distributions, the gas kinematics, and 
the quasar-host system dynamics with various models. Below are 
the main results: 


— The [C ir] and CO (9-8) lines and their underlying dust con- 
tinuum are all spatially resolved. The OH* emission is a 
point source, and the emission peak is only at a level of ^5c-. 
The line widths of the [C mn] and CO (9—8) emission are con- 
sistent, which suggests that these two lines trace the same 
gravitational potential. We observed asymmetric spectral line 
profiles for the [Cn] and CO (9-8) lines. The enhanced in- 
tensity of the lines toward the red part of the spectrum may 
be due to an asymmetric distribution of gas, or because the 
temperatures are higher on one side of the galaxy. 

— We used 2D elliptical Sérsic models to reproduce the ob- 
served intensity maps of the [Cu] and CO (9-8) lines and 
the dust continuum. The extended FIR dust continuum and 
the CO (9-8) emission can be fitted with a Sérsic distribu- 
tion with a Sérsic index of ~1 as well as an additional central 
point component. The point and extended components may 
represent emission from an AGN dusty molecular torus and 
the quasar host galaxy, respectively. The derived flux ratio 
from the AGN dust torus and the star-forming activity for 
the dust emission are consistent with our rest-frame UV-to- 
FIR SED decomposition result. The [C rr] emission follows a 
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flatter distribution (Sérsic index of 0.59), which may be due 
to high dust opacity increasing the FIR background of the 
[Cn] line and reducing the [C un] emission at smaller radii. 
The dust temperature drops with distance from the center. 
The effective radius of the dust continuum is smaller than 
that of the line emission and dust mass surface density, but 
is consistent with that of the SFR surface density. This may 
indicate that dust emission is a less robust tracer of the dust 
and gas distribution but is a decent tracer of the obscured 
star formation activity. The effective radius of the [C m] un- 
derlying dust continuum is consistent within ~10% with that 
of the CO (9-8) underlying dust continuum. The two wave- 
lengths are very close and shortward of the Rayleigh-Jeans 
tail; therefore, both are dominated by dust at similar temper- 
atures. The half-light radius of the [C n] line is smaller than 
that of the CO (9-8) line, which may be a reflection of the ra- 
dial dependence of the gas excitation and the radiation field. 
A larger [C m]-FIR deficit is observed in the center compared 
to the outer region, which is likely due to a higher dust tem- 
perature and higher dust opacity in the galaxy center. The 
spatially resolved KS relation (Zsgg ~ X7, ic) with Xy... 


measured from the CO (9-8) line of J231041855 i is steeper 
than in the local Universe and z ~ 2 — 4 starburst samples. 
However, when Xm, is measured from the [Cn] line, the 


relation (Xspr ~ DE 3740. 07) i 1s consistent with the local and 
low-redshift samples. 


A rotating gas disk is observed in both the [C ir] and CO (9- 
8) lines, as seen from the obvious gradients in the veloc- 
ity maps and the "S" shaped position-velocity diagrams. We 
applied 3D tilted ring models to both line data cubes with 
?DBAnoro. These two lines show similar kinematical signa- 
tures: the gas disk geometry (i.e., consistent inclination an- 
gle and position angle of the gas disk) and the gas motions 
(i.e., rotation dominated). The rotation velocities and veloc- 
ity dispersions traced by [C ir] and CO (9-8) lines are consis- 
tent. The circular velocities of [C n] and CO (9-8) lines are 
in excellent agreement, which is consistent with their iden- 
tical line width, indicating that they trace the same gravita- 
tional potential. We suggest that the [C n] and CO (9-8) gas 
are coplanar and corotating in this quasar host galaxy. The 
velocity dispersions and the V,../o ratio of our quasar host 
galaxy are comparable with other high-z populations, for ex- 
ample z > 6 quasar hosts and z ~ 2 star-forming galaxies. 
The OH* line shows a P-Cygni profile with an absorption at 
^—400 km/s, which may indicate an outflow with a neutral 
gas mass of (6.2 + 1.2) x 105 M, along the line of sight. 


For the purpose of quantifying the dynamical contributions 
from different matter components, we decomposed the cir- 
cular rotation curve measured from the high-resolution [C ri] 
line into four components (black hole, stars, gas, and dark 
matter) The whole quasar-host system appears baryonic 
matter dominated. The gas component is well constrained, 
which reveals a large quantity of gas with a mass on the order 
of 10!? M5. The stellar and dark matter components are not 
well constrained and depend on whether we allow the black 
hole mass to be a free parameter. We measured the black 
hole mass to be 2.971031 x 10° Mo. This is the first time that 
a SMBH mass has been measured dynamically at z ~ 6. This 
result is consistent with that measured from Mgu and C iv 
lines with the local scaling relations. The stellar component 
has a dynamical mass measured on the order of 10? Mo when 
the Universe was only ~0.93 Gyr old. This result is more ro- 
bust when we fix the black hole mass to be 1.8x 10? Mo from 
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the quasar spectrum. We note that we may miss the stellar 
components outside the gas and dust-emitting regions. The 
relations between the black hole mass and the stellar veloc- 
ity dispersion and the baryonic mass of this quasar indicate 
that the central SMBH evolved earlier than its host galaxy. 


In the future we plan to improve our understanding of 
the ISM distribution, gas kinematics, and system dynamics of 
quasars at the earliest epoch, with higher-resolution observations 
toward various rest-frame FIR ISM tracers with ALMA and to- 
ward the rest-frame optical and NIR stellar emission with the 
JWST. 
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Appendix A: The [C i] channel maps 


We present the [C n] and CO (9-8) channel maps of J2310+1855 in Figs. A.1 and A.2. The [C n] emission peak moves in a circular, 
counterclockwise path from 219 to —281 km/s, a signature of a rotating gas disk. 
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Fig. A.1. [Cu] channel maps of J2310+1855. The velocity takes the [C n] redshift from Wang et al. (2013) as a reference. The channel width is 
~ 17 km/s (corresponding to 15.625 MHz). The contour levels are [—2, 2, 4, 6, 8, 10, 12, 14] x rms, where rms = 0.17 mJy/beam. The white plus 
sign represents the HST quasar position. The [C rr] emission peak moves in a circular, counterclockwise path from 219 to —281 km/s. 
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Fig. A.2. CO (9-8) channel maps of J2310+1855. The velocity takes the [C n] redshift from Wang et al. (2013) as a reference. The channel width 
is ~ 32 km/s (corresponding to 15.625 MHz). The contour levels are [—2, 2, 4, 6, 8, 10, 12, 14] x rms, where rms = 0.10 mJy/beam. The white 
plus sign represents the HST quasar position. 


Appendix B: Two-dimensional elliptical Sérsic function 


The 1D Sérsic profile as a function of radius, r, is given by Sérsic (1963): 


X(r) = Xo expo", (B.D 
where X is the central surface brightness, r, is the scale length, which is the radius at which the surface brightness X(r) drops by 
el, and n is the shape parameter that controls the degree of curvature of the profile, the so-called Sérsic index. When n is large 
(i.e., > 4), rs is too small to measure; therefore, one often redefines the profile in terms of the half-light radius, re, also known as the 
effective radius: 

X(r) = Lp exp!) is (B.2) 


where the dependent variable k is coupled to n, and re = k”r;. k satisfies IT(2n) = B Yon k), where and y are the Gamma function 


: ; : 4 131 2194697 
and lower incomplete Gamma function, respectively. k can be approximated as 2n — i + 305 + EM - + ijasi754 + 39600717750, fOr 


n > 0.36 (Ciotti & Bertin 1999). 
The 2D elliptical Sérsic function can be expressed as 


f(xy) =A exp I4) +250 x0 )(y- yo) +e(-Yo) 1! ni (B.3) 
cos?0  sin?8 
a= UB + ud (B.4) 
x y 
sin20  sin20 
cx rue ee EEE A B.5 
b= OR R o 
: 2 2 
sin ð  cos^0 
c= dà + he” (B.6) 


Article number, page 29 of 33 


202302.00239v1 


chinaXiv 


A&A proofs: manuscript no. AA202244610 


where the coefficient A is the height of the peak, (xo, yo) is the center, hy, hy are the scale heights along the x and y axis before 
rotating by 6, and @ is the rotation angle in radians, which increases counterclockwise counting from the positive direction of the x 
axis, and n is the Sérsic index. We show surface plots, projected filled contour plots and surface brightness distributions in Fig. B.1 
for 2D elliptical Sérsic images with different rotation angles and Sérsic indices. 
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Fig. B.1. Surface plots, projected filled contour plots, and surface brightness distributions for 2D elliptical Sérsic images with (xo, yo) = (0, 0), 
h, = 4, and h, = 2. 0 and n are 7/3 and 0.5 (top left), 2/2 and 1 (top right), 27/3 and 2 (bottom left). The surface brightness as a function of radius 
shown in the bottom-right panel is measured using elliptical annuli with the same rotation angle of these images. The radius is the one along the 
major axis. The black, red and blue dots are measured from the top-left, top-right, and bottom-left images, respectively. The black, red, and blue 
lines are the fits to these surface brightness distributions with 1D Sérsic function in Eq. B.1. The best-fitted Sérsic indices are labeled in the same 
colors as the fitting lines. 


We use a 2D elliptical Sérsic model to fit the CO (9-8) line and the [C ir] and CO (9-8) underlying dust continuum. And as 
motivated by the existence of AGN torus, we also fit these intensity maps with an additional unresolved nuclear component to 
represent the dusty and molecular torus. The results are shown in Figs. 7, B.2, and B.3. 


Appendix C: Tilted ring model 


In this section, we illustrate the parameters used in the tilted ring model. We only observe the line-of-sight velocities as a function of 
position on the sky, which cannot be directly used to investigate the mass distributions in galaxies. We need to measure the rotation 
velocity Vot as a function of distance R to the center (the rotation curve). If we assume that the gas is rotating symmetrically in a 
disk, a number of tilted ring model parameters can define the observed velocity field of a galaxy. The tilted ring model (see Fig. C.1) 
consists of a set of concentric rings characterized by various geometrical and kinematical components. The geometrical components 
are: 

* Xo, Yo: Sky coordinates of the rotation center of the galaxy. 

e i (R): the inclination (i.e., the angle between the normal to the plane of the galaxy and the line-of-sight. 

e ¢ (R): the position angle of the major axis of a ring projected onto the sky (i.e., an ellipse). This is an angle taken in a 
counterclockwise direction from the north of the sky to the major axis of the receding half of the galaxy. 
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Fig. B.2. As in Fig. 7 but for the [C ir] underlying dust continuum. The peak values and rms within the dashed gray square are 0.101 and 0.028, 
and 0.099 and 0.028 Jy/beam km/s for the Sérsic and P+Sérsic residual maps, respectively. The shape of the synthesized beam with a FWHM size 
of 07113 x 07092 is plotted as white ellipses. We measured the [Cm] underlying dust continuum luminosity surface density using elliptical rings 
with the ring width along the major axis of half (0705) that of the major axis of the [Cu] clean beam size, the rotation angle equal to PA (=199°) 
and the ratio of semiminor and semimajor axis — b/a of cos(i) (i = 42°), where PA and i come from the [Cn] line kinematic modeling (listed in 
Table 3). 


e ¢ (R): the azimuthal angle measured between the major axis on the sky plane and the line connecting the galaxy center (xo, yo) 
and the data point (X, Y) in the galaxy plane. It is related to i (R), ø (R), (xo, yo). 


The kinematical components are: 
e Vys: the velocity of the center of the galaxy with respect to the Sun, the so-called systemic velocity. 
* Vior(R): the rotation velocity at distance R from the center. This velocity is assumed to be constant in a ring. 


The V;ot is perpendicular to the ring and has components in the X and Y directions. Vx (= Vrot sin Z) is parallel to the xy plane and 
has no component in the z direction. However, Vy (= Vis: cos Z) has a component in the z direction: V; = Vy sin i. This is the radial 
velocity (Vios) that we observe at a point (x, y). It must be corrected for the systemic (line-of-sight) velocity V... Thus, 


Vios(x, y) = Vsys + Vz = Vsys + Vrot cosZ sini (C.1) 
ere = me PI ese (C.2) 
R R 
. , Y  -(Q-xcosó - (y — yo) sing 
sing = p Roosi : (C.3) 


Here, we also put the PBARoLo modeling result of CO (9-8) line in Fig. C.2. 
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Fig. B.3. As in Fig. 7 but for the CO (9-8) underlying dust continuum. The peak values and rms within the dashed gray square are 0.043 and 0.014, 
and 0.041 and 0.014 Jy/beam km/s for the Sérsic and P+Sérsic residual maps, respectively. The shape of the synthesized beam with a FWHM size 
of 07192 x 0/156 is plotted as a white ellipse in each panel. 


X (East) © X (major axis) 


Fig. C.1. Simple schematic diagram for the tilted ring model. The pink semiellipse represents the galaxy plane with a coordinate system of XYZ. 
The red-based color lines are on or perpendicular to the galaxy plane. The light-blue ellipse presents the sky plane with a coordinate system of 
xyz. The blue color lines are on or perpendicular to the sky plane. We takes the x direction as the east and the y direction as the north. The major 
axis on the sky plane coincides with the X axis on the galaxy plane. The line-of-sight direction is along the z direction from the bottom to the 
top. The galaxy plane is inclined with an angle of i (the inclination angle, i.e., the angle between the line of sight and the perpendicular line to the 
galaxy plane XY). The data point (X, Y) with a distance of R to the center rotates counterclockwise with a rotation velocity of V;4, on the galaxy 
plane. What we can see is the line-of-sight component V, of the rotation velocity Vo. The angle between the major axis on the sky plane and the 
line connecting the galaxy center (xo, yo) and the data point (X, Y) is the azimuth angle Z. The position angle $ is defined as the angle counting 
from the north (y axis) to the major axis of the receding half on the galaxy plane. After tilting, the data point (X, Y), the distance R and the azimuth 
angle Z on the galaxy plane project to (x, y), R’ and £ on the sky plane. 
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Fig. C.2. Similar to that of Fig. 10, but for the CO (9-8) line kinematic modeling. Upper panels — Shape of the synthesized beam with a FWHM 
size of 07187 x 07153 is plotted in the bottom-left corner of each panel. Lower panels — Left and middle panels: Contour levels are [—2, 2, 4, 8, 
16] x 0.093 mJy/beam for both the data and the model. Right panel: The spectral resolution is 31.25 MHz, corresponding to 64 km/s. 
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